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PREFACE 


1.  The  Corps  of  Engineers,  through  its  Civil  Works  program,  has  spon- 
sored, over  the  past  23  years,  research  into  the  behavior  and  characteris- 
tics of  tidal  inlets.  The  Corps'  interest  in  tidal  inlet  research  stems 
from  its  responsibilities  for  navigation,  beach  erosion  prevention  and 
control,  and  flood  control.  Tasked  with  the  creation  and  maintenance  of 
navigable  U.S.  waterways,  the  Corps  dredges  millions  of  cubic  yards  of 
material  each  year  from  tidal  inlets  that  connect  the  ocean  with  bays, 
estuaries,  and  lagoons.  Design  and  construction  of  navigation  improve- 
ments to  existing  tidal  inlets  are  an  important  part  of  the  work  of  many 
Corps'  offices.  In  some  cases,  design  and  construction  of  new  inlets  are 
required.  Development  of  information  concerning  the  hydraulic  character- 
istics of  inlets  is  important  not  only  for  navigation  and  inlet  stability, 
but  also  because  inlets,  by  allowing  for  the  ingress  of  storm  surges  and 
egress  of  flood  waters,  play  an  important  role  in  the  flushing  of  bays 
and  lagoons. 

2.  A research  program,  the  General  Investigation  of  Tidal  Inlets 
(GITI),  was  developed  to  provide  quantitative  data  for  use  in  design 

of  inlets  and  inlet  improvements.  It  is  designed  to  meet  the  following 
objectives : 

To  determine  the  effects  of  wave  action,  tidal  flow,  and  related 
phenomena  on  inlet  stability  and  on  the  hydraulic,  geometric,  and 
sedimentary  characteristics  of  tidal  inlets;  to  develop  the  knowl- 
edge necessary  to  design  effective  navigation  improvements,  new 
inlets,  and  sand  transfer  systems  at  existing  tidal  inlets;  to 
evaluate  the  water  transfer  and  flushing  capability  of  tidal 
inlets;  and  to  define  the  processes  controlling  inlet  stability. 

3.  The  GITI  is  divided  into  three  major  study  areas:  (a)  inlet 
classification,  (b)  inlet  hydraulics,  and  (c)  inlet  dynamics. 

a.  Inlet  Classification.  The  objectives  of  the  inlet  classifi- 
cation study  are  to  classify  inlets  according  to  their  geometry,  hydrau- 
lics, and  stability,  and  to  determine  the  relationships  that  exist  among 
the  geometric  and  dynamic  characteristics  and  the  environmental  factors 
that  control  these  characteristics.  The  classification  study  keeps  the 
general  investigation  closely  related  to  real  inlets  and  produces  an 
important  inlet  data  base  useful  in  documenting  the  characteristics  of 
inlets . 


b.  Inlet  Hydraulics.  The  objectives  of  the  inlet  hydraulics 
study  are  to  define  tide-generated  flow  regime  and  water  level  fluctua- 
tions in  the  vicinity  of  coastal  inlets  and  to  develop  techniques  for 
predicting  these  phenomena.  The  inlet  hydraulics  study  is  divided  into 
three  areas:  (1)  idealized  inlet  model  study,  (2)  evaluation  of  state- 
of-the-art  physical  and  numerical  models,  and  (3)  prototype  inlet 
hydraulics . 
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(1)  The  Idealized  Inlet  Model.  The  objectives  of  this 
model  study  are  to  determine  the  effect  of  inlet  configurations  and 
structures  on  discharge,  head  loss  and  velocity  distribution  for  a num- 
ber of  realistic  inlet  shapes  and  tide  conditions.  An  initial  set  of 
tests  in  a trapezoidal  inlet  was  conducted  between  1967  and  1970.  How- 
ever, in  order  that  subsequent  inlet  models  are  more  representative  of 
real  inlets,  a number  of  "idealized"  models  representing  various  inlet 
morphological  classes  are  being  developed  and  tested.  The  effects  of 
jetties  and  wave  action  on  the  hydraulics  are  included  in  the  study. 

(2)  Evaluation  of  State-of-the-Art  Modeling  Techniques. 

The  objectives  of  this  part  of  the  inlet  hydraulics  study  are  to  deter- 
mine the  usefulness  and  reliability  of  existing  physical  and  numerical 
modeling  techniques  in  predicting  the  hydraulic  characteristics  of  inlet- 
bay  systems,  and  to  determine  whether  simple  tests,  performed  rapidly  and 
economically,  are  useful  in  the  evaluation  of  proposed  inlet  improvements. 
Masonboro  Inlet,  North  Carolina,  was  selected  as  the  prototype  inlet  which 
would  be  used  along  with  hydraulic  and  numerical  models  in  the  evaluation 
of  existing  techniques.  In  September  1969  a complete  set  of  hydraulic 
and  bathymetric  data  was  collected  at  Masonboro  Inlet.  Construction  of 
the  fixed-bed  physical  model  was  initiated  in  1969,  and  extensive  tests 
have  been  performed  since  then.  In  addition,  three  existing  numerical 
models  were  applied  to  predict  the  inlet's  hydraulics.  Extensive  field 
data  were  collected  at  Masonboro  Inlet  in  August  1974  for  use  in  evaluat- 
ing the  capabilities  of  the  physical  and  numerical  models. 

(3)  Prototype  Inlet  Hydraulics.  Field  studies  at  a number 
of  inlets  are  providing  information  on  prototype  inlet-bay  tidal  hydraulic 
relationships  and  the  effects  of  friction,  waves,  tides,  and  inlet  mor- 
phology on  these  relationships. 

a.  Inlet  Dynamics.  The  basic  objective  of  the  inlet  dynamics 
study  is  to  investigate  the  interactions  of  tidal  flow,  inlet  configura- 
tion, and  wave  action  at  tidal  inlets  as  a guide  to  improvement  of  inlet 
channels  and  nearby  shore  protection  works.  The  study  is  subdivided 
into  four  specific  areas:  (1)  model  materials  evaluation,  (2)  movable- 
bed  modeling  evaluation,  (3)  reanalysis  of  a previous  inlet  model  study, 
and  (4)  prototype  inlet  studies. 

(1)  Model  Materials  Evaluation.  This  evaluation  was  initi- 
ated in  1969  to  provide  data  on  the  response  of  movable-bed  model  mate- 
rials to  waves  and  flow  to  allow  selection  of  the  optimum  bed  materials 
for  inlet  models. 

(2)  Movable-Bed  Model  Evaluation.  The  objective  of  this 
study  is  to  evaluate  the  state-of-the-art  of  modeling  techniques,  in 
this  case  movable-bed  inlet  modeling.  Since,  in  many  cases,  movable-bed 
modeling  is  the  only  tool  available  for  predicting  the  response  of  an 
inlet  to  improvements,  the  capabilities  and  limitations  of  these  models 
must  be  established. 
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(3)  Reanalysis  of  an  Earlier  Inlet  Model  Study.  In  1955, 
a report  entitled,  "Preliminary  Report:  Laboratory  Study  of  the  Effect 
of  an  Uncontrolled  Inlet  on  the  Adjacent  Beach,"  was  published  by  the 
Beach  Erosion  Board  (now  CERC) . A reanalysis  of  the  original  data  is 
being  performed  to  aid  in  planning  of  additional  GITI  efforts. 

(4)  Prototype  Dynamics.  Field  and  office  studies  of  a 
number  of  inlets  are  providing  information  on  the  effects  of  physical 
forces  and  artificial  improvements  on  inlet  morphology.  Of  particular 
importance  are  studies  to  define  the  mechanisms  of  natural  sand  bypassing 
at  inlets,  the  response  of  inlet  navigation  channels  to  dredging  and 
natural  forces,  and  the  effects  of  inlets  on  adjacent  beaches. 

4.  This  report  was  directed  toward  the  evaluation  of  the  state-of- 
the-art  of  numerical  and  physical  models  and  modeling  techniques  for 
tidal  inlets.  However,  technology  advanced  so  rapidly  during  the  study 
period  and  the  review  process  that  a current  state-of-the-art  report  did 
not  result.  Nevertheless,  this  report  is  believed  a significant  contri- 
bution toward  the  evaluation  of  current  modeling  techniques  and  will  be 
useful  to  coastal  engineers  for  several  years.  The  evaluation  was  con- 
ducted by  calibrating  a physical  model  and  three  numerical  models  with 
data  from  the  prototype.  Since  no  verification  data  were  available  from 
the  prototype,  the  models  were  not  verified.  However,  the  models  were 
all  applied  to  significantly  different  inlet  hydrographic  conditions 
and  the  results  of  the  numerical  models  are  compared  with  those  of  the 
physical  model. 

5.  As  an  aid  in  model  comparison  and  evaluation,  the  two-dimensional 
equations  of  motion  which  govern  tidal  flows  in  inlets  (and  many  similar 
phenomena)  are  systematically  derived  from  the  Navier-Stokes  equations. 
This  mode  of  derivation  provides  a clear  indication  of  the  processes 
often  neglected  in  the  study  of  tidal  flows. 

6.  Reports  on  each  of  the  four  models  involved  are  published  as 
separate  appendixes  to  this  report.  A later  report  will  present  a 
comparison  of  new  prototype  observations  with  predictions  of  the  physical 
model,  and  of  the  best  numerical  model  for  prototype  conditions  not 

used  in  model  calibration. 
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CONVERSION  FACTORS,  U.S.  CUSTOMARY  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


U.S.  customary  units  of  measurement 
to  metric  (SI)  units  as  follows: 

used  in  this  report  can  be  converted 

Multiply 

by 

To  obtain 

inches 

25.4 

millimeters 

2.54 

centimeters 

square  inches 

6.452 

square  centimeters 

cubic  inches 

16.  39 

cubic  centimeters 

feet 

30.39 

centimeters 

0.3048 

meters 

square  feet 

0.0929 

square  meters 

cubic  feet 

0.0283 

cubic  meters 

yards 

0.9144 

meters 

square  yards 

0.836 

square  meters 

cubic  yards 

0.7646 

cubic  meters 

miles 

1.6093 

kilometers 

square  miles 

259.0 

hectares 

knots 

1.8532 

kilometers  per  hour 

acres 

0.4047 

hectares 

foot-pounds 

1.3558 

newton  meters 

millibars 

1.0197  x 10 

"3  kilograms  per  square  centimeter 

ounces 

28.35 

grams 

pounds 

453.6 

grams 

0.4536 

kilograms 

ton,  long 

1.0160 

metric  tons 

ton,  short 

0.9072 

metric  tons 

degrees  (angle) 

0.1745 

radians 

Fahrenheit  degrees 

5/9 

Celsius  degrees  or  Kelvins1 

1To  obtain  Celsius 

(C)  temperature 

readings  from  Fahrenheit  (F)  readings, 

use  formula:  C = 

(5/9)  (F  -32). 

To  obtain  Kelvin 

(K)  readings,  use 

formula:  K = (5/9)  (F  -32)  + 273.15. 

SYMBOLS  AND  DEFINITIONS 


amplitude  of  a sinusoidal  wave 
arbitrary  grid  point 


A(x,t)  = / D dy,  the  cross  section  of  an  inlet  channel 


(Sec.  V) 


cross-sectional  area  of  an  inlet  at  its  junction  with  a bay 
or  lagoon 

area  of  the  bay  connected  to  the  ocean  by  an  inlet 
Ac(x,t) , cross-sectional  area  of  an  inlet 


amplitude  of  the  n'th  sinusoidal  term  in  a Fourier  expansion 


cross  section  of  the  n'th  inlet  subchannel  (Sec.  V) 


cross  section  of  an  inlet  at  its  junction  with  the  sea 
(Sec.  V) 

turbulent  mixing  coefficients  in  a horizontal  plane 


1 

g(Rg  - fig,) 


/ 


*s 


an  area  characterizing  an  inlet 


coefficients  of  the  quadratic  equation  in  standard  form 
(Sec.  V) 

coefficient  of  the  cosine  term  of  the  n'th  harmonic  in  a 
Fourier  expansion 

a second  arbitrary  grid  point 

as  a subscript,  indicates  that  the  primary  variable  pertains 
to  the  bay 

amplitude  of  the  n'th  sinusoidal  term  in  a second  Fourier 
expansion 

b(x,t)  = Y2  ~ Yl  > width  of  inlet  at  water  surface 
width  of  the  m'th  subchannel  in  an  inlet 


SYMBOLS  AND  DEFINITIONS— Continued 


coefficient  of  the  sine  term  of  the  n'th  harmonic  in  a 
Fourier  expansion 

c phase  velocity  of  a wave 

Cnax  maximum  phase  velocity  of  a wave 

D D(x,y,z,t),  total  depth  of  fluid  layer 

Dm  D applied  to  the  model  (Sec.  Ill) 

depth  of  the  m'th  subchannel  in  an  inlet  (Sec.  V) 

Umax  maximum  value  of  the  total  depth 
D min  minimum  value  of  the  depth 

Dp  D applied  to  prototype 

D cross-section  average  value  of  D 

D*  dimensionless  form  of  D 

D'  D-D 

F arbitraiy  function 

f 2ft  sin  <(>,  Coriolis  parameter 

f2  2ft  cos  <f> , second  Coriolis  parameter 

g acceleration  of  gravity 

height  of  the  barometer  column 

h h(x,y,t),  elevation  of  top  of  a fluid  layer,  usually  the  water 

surface 

h as  a subscript,  indicates  that  the  variable  is  evaluated  at  the 

free  surface  of  the  water 

hp  h applied  to  the  prototype 

hj  mean  elevation  of  the  water  surface  over  a period  much  longer  than 

that  used  in  defining  h;  hj  is  generally  interpreted  as  mean  sea 
level 

h mean  value  of  h over  a few  wave  periods 

h as  a subscript,  indicates  that  the  function  is  evaluated  at  the 

mean  position  of  the  free  surface 
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SYMBOLS  AND  DEFINITIONS--Continued 


h' 

h' 


h 

B>  ^5 

h* 

h 


j 

K 

Ke 

KHW 

K m 

h 

L 

LJ 

M 


m 

m 

N 

n 


h - h 

h'(x,y,t)  = h - ft,  deviation  of  water  surface  elevation  from 
the  cross-channel  average 

= , y 2 - 

h(x,t)  = b-1  I h dy,  average  value  of  h across  an  inlet 

y\ 

channel 

values  of  ii  at  bay  and  seaward  ends  of  an  inlet 
dimensionless  form  of  h 

the  contribution  of  waves  to  the  elevation  of  the  water  surface 

/T 

arbitrary  integer 

dimensionless  friction  coefficient 
Keulegan's  coefficient  of  repletion 

Keulegan's  (1967)  coefficient  of  repletion  as  defined  by 
Huval  and  Wintergerst  (1977) 

the  value  of  K that  applies  to  the  model 

friction  coefficient  of  the  n'th  subchannel  in  an  inlet 

value  of  K that  applies  to  the  prototype 

length  of  a gravity  wave  or  other  disturbance 

length  of  inlet 

maximum  value  of  m,  when  m is  used  as  a running  index 

coefficient  of  inertia  term  in  Keulegan's  inlet  equations 
(Sec.  V) 

a subscript  indicating  a running  index  (Sec.  V) 

a subscript  to  indicate  that  the  primary  variable  applies  to  a 
model  (Sec.  Ill) 

maximum  value  of  n in  summation 
distance  along  a wave  crest  (Sec.  Ill) 
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SYMBOLS  AND  DEFINITIONS--Continued 


n Manning's  n,  a dimensionless  resistance  coefficient 

n running  index  in  summations 

n as  a subscript  indicates  that  the  primary  variable  pertains  to 

the  prototype 

n'  second  running  index 

p p(x,y,z,t),  pressure 

p^  pressure  at  the  top  of  a fluid  layer 
p pressure  at  the  bottom  of  a fluid  layer 

p average  value  of  p over  a small  volume  on  a time  duration  of 

a few  wave  periods 

pn  time-mean  pressure  at  the  free  surface  of  the  fluid 

p'  contribution  of  the  perturbation  pressure  to  p 

p contribution  of  wave  motion  to  p 

Q uA,  inlet  discharge  in  a one-dimensional  flow 

Q (U2  + V2)1//2,  transport  magnitude  in  a vertical  column  having 

unit  cross-sectional  area  in  two-dimensional  flow 


q 

^ in 
% 


discharge  in  the  n'th  subchannel 

horizontal  velocity  in  the  direction  of  s 

flow  into  bay  which  does  not  pass  through  inlet  (Sec.  V) 

time-mean  velocity  at  the  outer  edge  of  the  viscous  sublayer 
(Sec.  II) 


time-mean  vertically  averaged  current  velocity  (Sec.  II) 
q contribution  of  waves  to  q 

qz  q at  z = Z 

R rainfall  rate;  specifically,  the  net  excess  of  rainfall,  runoff 

and  bottom  seepage  into  the  basin  over  evaporation,  and  bottom 
seepage  out  of  the  basin 

RMSP  standard  deviation  of  prototype  measurements  (eq.  112) 
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SYMBOLS  AND  DEFINITIONS--Continued 


RMSA  root-mean-square  difference  between  two  measurements  (eq.  Ill) 

Re  closure  term  for  the  one-dimensional  channel  equations  defined 

by  equation  (74)  in  Section  V 

Rg  v/f,  radius  of  the  inertial  circle 

S as  a subscript,  indicates  that  the  primary  variable  is  to  be 

evaluated  at  x = Xg,  the  inlet  end  nearest  the  sea 

Sxo>  Syo  closure  terms  for  the  two-dimensional  momentum  equations, 
other  than  S xt,  $yt>  $xv>  Syy»  S xu,  or  Syu 

^xt>  Syt  closure  terms  representing  molecular  and  turbulent  stresses 
in  two-dimensional  momentum  equations 

S^ry,  ^yv  closure  terms  arising  from  vertical  integration  for  the 
two-dimensional  momentum  equations 

S xu,  S yW  closure  terms  arising  from  surface  waves  for  the  two- 

dimensional  momentum  equations 

Srr,  S^,  components  of  the  Reynolds  stress  tensor  due  to  waves 

S yx>  $yy  (eqs.  16,  17,  and  18) 

s arbitrary  horizontal  direction 

s distance  in  the  direction  of  wave  propagation 

s.  slip  coefficient  in  the  Jelesnianski  (1967)  bottom  friction 

law 

T period  of  gravity  wave;  may  be  surface  wave  or  tide  wave 

t time 

tm  t applied  to  the  model 

t p t applied  to  the  prototype 

t*  dimensionless  form  of  t 

h 

U U(x,y,t)  = f udz,  transport  parallel  to  the  x-axis  in  a 

Z 

vertical  column  having  unit  cross-sectional  area 
U applied  to  the  model 
Up  U applied  to  the  prototype 
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SYMBOLS  AND  DEFINITIONS--Continued 


U* 


u,  v,  w 


U,  V,  w 

u'  , v' , w' 
U,  V,  w 

uh>  vh 
uh>  vh 
“z>  VZ 


UB>  US 


°n 

u" 

u",  v" 


dimensionless  form  of  U 

components  of  the  total  velocity  parallel  to  the  x- , y- , 
and  z-axes 

elementary  averages  of  the  velocity  components  parallel  to 
the  x-,  y- , and  z-axes 

part  of  the  velocity  components  due  to  turbulence 
part  of  the  velocity  components  due  to  surface  waves 
horizontal  velocity  components  evaluated  at  z = h 

horizontal  velocity  components  evaluated  at  z = h 

horizontal  velocity  components  evaluated  at  z = Z 

. y 2 

u(x,t)  = A^1  J U dy,  mean  value  of  u in  an  inlet  cross 
Vl 

section 

mean  horizontal  velocity  along  the  inlet  axis  at  bay  and 
seaward  ends  of  inlet 

mean  velocity  in  the  n'th  subchannel  of  an  inlet 
u"(x,y,z,t)  = u(x,y,z,t)  - 5(x,t) 

departures  of  the  horizontal  velocity  components  from  the 
vertically  average  values 


V 


h 

V(x,y,t)  = / vdz,  transport  parallel  to  the  y-axis  in  a 
Z 

vertical  column  having  unit  cross-sectional  area 
volume  of  bay 

} ?B  Y 2 

V|  - J J j (R  + E)  dy  dx  (eq.  85) 

9 XS  Yi 

V applied  to  the  model 

V applied  to  the  prototype 
dimensionless  form  of  V 
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SYMBOLS  AND  DEFINITIONS--Continued 


Vf-  volume  of  inlet 

v horizontal  velocity  in  arbitrary  direction 

Mm  wind  velocity  for  the  model 

Vp  wind  velocity  for  the  prototype 

Y 2 

v v(x,t)  = Ac1  / V dy,  mean  value  of  v in  an  inlet  cross 

yi 

section,  negligibly  small  in  most  cases 
w^  vertical  velocity  at  top  of  fluid  layer  (App.  B,  eq.  B15) 

w7  vertical  velocity  at  bottom  of  fluid  layer  (App.  B,  eq.  B14) 

w2  vertical  velocity  due  to  waves  at  an  arbitrary  value  of  z 

wJ^  vertical  velocity  at  the  free  surface  as  a result  of 

turbulence 

X arbitrary  statistical  function 

x,  y,  z principal  axes  in  a right-handed  coordinate  system;  z is 

positive  upward 

xB,  Xg  value  of  x at  the  bay  (ocean)  end  of  an  inlet 

Xm,  ym,  zm  x,  y,  z applied  to  model 

Xp»  Yp>  zp  x,  y,  z applied  to  prototype 

x*,  y* , z*  dimensionless  forms  of  x,  y,  z 

Y second  arbitrary  statistical  function 

yj,  y2  values  of  y at  opposite  ends  of  an  inlet  cross  section 

Z Z(x,y,t),  elevation  of  the  bottom  of  a fluid  layer 

Z as  a subscript,  Z indicates  that  the  primary  variable 

should  be  evaluated  at  z = Z 

Zj  arbitrary  value  of  z < hj 
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SYMBOLS  AND  DEFINITIONS--Continued 


a a^/oip,  horizontal  scale  ratio 

am  dimensional  factor  for  the  horizontal  dimension  in  a model 

otp  dimensional  factor  for  the  horizontal  dimension  in  the  prototype 

cij  arbitrary  exponent 

6 ^/Tj/Sp , vertical  scale  ratio 

Bm  dimensional  factor  for  the  vertical  dimension  in  the  model 

Bp  dimensional  factor  for  the  vertical  dimension  in  the  prototype 

Y y = +1  indicates  floodtide;  y = -1  indicates  ebbtide  (Sec.  V) 

Yffl  dimensional  factor  for  time  in  the  model 

Y dimensional  factor  for  time  in  the  prototype 

r 

Ah  increment  in  height 

As  increment  of  horizontal  distance  in  any  direction 

At  increment  of  time 

Ax  increment  of  distance  along  x-axis 

6 thickness  of  the  viscous  sublayer 

n a value  of  z,  h < n < h,  used  in  an  application  of  the  mean  value 
theorem  of  integral  calculus 

9 arbitrary  angle 

9 direction  of  a wave  ray 

9n  phase  angle  in  a Fourier  expansion 

^ horizontal  length  of  an  underwater  slope,  used  only  in  the 
discussion  of  physical  models 

y coefficient  of  dynamic  molecular  viscosity 

v y/P,  coefficient  for  kinematic  viscosity 

v effective  viscosity  coefficient  in  the  Jelesnianski  (1967)  friction 
model 
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SYMBOLS  AND  DEFINITIONS--Continued 
3.14159  


P 


fluid  density 


p"  p"  = P - P 

h 

p = (h  - D)”1  j p dz 
Z 


xx> 

xyx> 

xzx> 

\xy> 

Tyy> 

Tzy> 

XZ’ 

xyz> 

Tzz 

components  of  the  stress  tensor.  The  first  subscript 
indicates  the  face  of  an  elementary  cube  on  which  the 
stress  acts;  the  second  subscript  indicates  the 
direction  in  which  the  stress  acts 


(xzx)z>  (xzy^Z  stresses  in  a h^  at  the  bottom  of  a fluid  layer 
parallel  to  the  x-  and  y-axes 


(Tzx)h>  ^zy^h  stresses  in  a horizontal  plane  at  the  upper  boundary 
of  a fluid  layer  parallel  to  the  x-  and  y-axes 

(.tzx) Z mean  bottom  stress  at  an  inlet  cross  section 

<t>  latitude 


Y 

Y 
? 
$ 

SI 


u 


phase  angle  for  a second  Fourier  expression 
angle  between  the  x-axis  and  the  east  direction 
arbitrary  function 


V - ip 


value  of  ip  averaged  over  a small  interval  of  time 
and  space 

angular  speed  of  the  earth's  rotation  = 7.292  * 10-5 
radians  per  second 

w(x,y),  weighting  function  for  velocity  along  the  axes 
of  an  inlet  in  lumped  parameter  model  = u/Q^  = A^/A^(x) 
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COMPARISON  OF  NUMERICAL  AND  PHYSICAL  HYDRAULIC  MODELS, 
MASONBORO  INLET,  NORTH  CAROLINA 


by 

D.  Lee  Harris 
and 

B.R.  Bodine 
I.  INTRODUCTION 


1 . Historical  Perspective. 

In  1969,  the  Coastal  Engineering  Research  Board  (CERB)  recommended 
that  a comparative  test  of  the  accuracy  of  numerical  models  be  conducted 
in  conjunction  with  a distorted  scale,  fixed-bed  model  already  planned 
for  Masonboro  Inlet,  North  Carolina.  The  planned  fixed-bed  physical 
model  was  to  be  used  to  predict  changes  in  inlet  hydraulic  characteris- 
tics resulting  from  improvements  at  the  inlet  entrance.  A bathymetric 
survey  of  the  inlet  was  made  in  September  1969.  Tide  and  current  obser- 
vations through  one  tidal  period  were  collected  on  12  September  1969  at 
the  locations  shown  in  Figure  1 to  provide  information  needed  for  cali- 
brating the  fixed-bed  model. 

Implementation  of  the  CERB  recommendation  was  initiated  in  1971  with 
the  issuance  of  a Request  for  Proposals  (bids)  (App.  A).  In  addition, 
the  bidders  were  supplied  with  copies  of  the  data  collected  in  1969  for 
use  in  calibrating  the  numerical  models.  The  plan  called  for  model  cali- 
bration with  the  1969  hydrographic  and  hydraulic  data.  All  of  the  models 
were  to  be  applied  to  November  1964  and  June  1967  hydrographic  data, 
using  standard  ocean  tides  to  drive  the  models.  The  results  of  the 
application  of  the  models  to  the  new  conditions  would  then  be  compared 
to  evaluate  the  existing  capability  to  use  numerical  hydraulic  models  to 
predict  hydraulic  conditions  at  coastal  inlets. 

After  the  bids  were  received  and  the  work  under  resulting  contracts 
was  completed,  it  was  decided  that  the  physical  model  would  be  operated 
for  conditions  of  June  1966  rather  than  conditions  of  June  1967.  Since 
there  were  great  differences  in  the  hydrography  for  the  two  conditions, 
the  hydraulic  comparisons  would  have  little  meaning.  Therefore,  the 
comparisons  of  the  numerical  models  with  the  results  of  the  physical 
model  for  the  1967  conditions  are  not  discussed  in  this  report.  However, 
results  from  the  1967  hydrography  are  included  in  the  numerical  modeling 
reports  published  as  separate  appendixes  to  this  report:  Chen  and  Hembree 
(1977),  Huval  and  Wintergerst  (1977),  Masch,  Brandes , and  Reagan  (1977), 
and  Sager  and  Seabergh  (1977).  New  prototype  measurements  made  at  Mason- 
boro in  1974  were  compared  with  results  of  the  physical  model  study  by 
Seabergh  and  Mason  (1975).  This  comparison  of  the  1974  prototype  data 
and  numerical  model  results  will  be  published  in  a later  report. 

As  a result  of  the  Request  for  Proposals,  two  two-dimensional  models, 
constructed  along  the  lines  of  Reid  and  Bodine  (1968),  and  a lumped 
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Locations  of  tide  and  current  observations,  Masonboro  Inlet. 


parameter  model  were  selected  from  the  many  models  proposed  for  study. 

The  two-dimensional  numerical  models  were  expected  to  provide  tidal 
heights  and  mean  current  estimates  for  one  tidal  cycle  at  each  measure- 
ment location  (Fig.  1) . The  lumped  parameter  model  would  provide  only 
one  estimated  current  hydrograph  (a  cross-sectional  average  at  the  throat 
of  the  inlet)  and  one  estimated  tide  hydrograph,  indicative  of  average 
conditions  inside  the  inlet. 

The  use  of  both  physical  and  numerical  hydraulic  models  involves  a 
preliminary  stage  or  stages  during  which  the  relationships  between  the 
model  and  the  prototype  are  investigated  and  established.  Various  terms 
have  been  used  to  describe  this  process;  e.g.,  confirmation , adjustment, 
verification,  and  calibration.  These  terms  are  sometimes  used  inter- 
changeably in  the  hydraulic  model  literature  with  definitions  varying  by 
place  and  time.  The  following  definitions  are  used  in  this  study: 

(a)  Calibration.  The  process  where  a hydraulic  model  is  checked 
with  prototype  data  and  systematically  adjusted  to  reproduce  the  water 
level  and  current  data  from  corresponding  prototype  driving  forces. 

(b)  Verification.  The  process  where  independent  prototype  data 
(data  not  used  in  the  calibration)  are  used  to  verify  that  a calibrated 
hydraulic  model  produces  satisfactory  prototype  results. 

Since  data  from  the  prototype  are  available  for  only  one  condition 
(the  1969  hydrography  and  the  September  1969  tidal  height  and  tidal 
current  measurements),  the  models  used  in  this  study  could  only  be  cali- 
brated. It  was  assumed  at  the  start  of  this  study  that  more  confidence 
could  be  placed  in  the  physical  model  results  than  in  the  numerical  model 
results.  Therefore,  a comparison  of  the  numerical  model  results  with  the 
1964  results  of  the  physical  model  could  be  considered  as  verification  of 
the  numerical  models. 

The  concept  of  verification  of  numerical  models  by  comparison  with 
physical  models  was  abandoned  when  it  was  determined  that  one  of  the 
numerical  models  performed  as  well  as  the  physical  model  in  the  calibra- 
tion process  and  that  physical  reality  is  approximated  differently  by 
numerical  models  than  by  physical  models. 

2.  Types  of  Models. 

A model  may  be  a (a)  geometrically  similar  physical  replica  of  the 
prototype;  (b)  distorted  physical  replica  of  the  prototype  for  which 
appropriate  corrections  can  be  made;  (c)  set  of  mathematical  equations 
that  express  the  laws  which  govern  the  phenomena  of  interest  in  the  pro- 
totype; or  (d)  physically  dissimilar  system  governed  by  a similar  set  of 
equations,  together  with  a statement  of  the  relationships  between  the 
two  systems.  This  study  primarily  discusses  the  consideration  of  dis- 
torted physical  replicas  of  tidal  inlet  flow  and  a set  of  mathematical 
equations  which  describe  flow  within  and  near  the  inlet.  However,  un- 
distorted models  and  electrical  analogs  of  the  flow  could  be  used. 
Ishiguro  (1972)  reviewed  the  use  of  electrical  analogs  in  the  solution 
of  tidal  hydraulics  problems. 


23 


Although  physical  and  numerical  models  are  different  approximations 
of  the  same  prototype  phenomena,  it  is  possible  that  each  technique  is 
superior  to  the  other  in  some  respects.  One  model  might  be  more  useful 
for  investigating  specific  problems;  both  models  might  be  equally  satis- 
factory for  other  problems,  even  though  each  yields  slightly  different 
results . 

Although  laws  which  govern  a phenomenon  have  not  been  firmly  estab- 
lished, the  design  of  useful  physical  model  experiments  is  often  possible, 
provided  certain  dimensionless  parameters,  formed  from  the  most  important 
variables  associated  with  the  phenomenon,  have  the  same  value  for  both 
model  and  prototype.  This  procedure  is  widely  used  in  hydraulics,  and 
is  discussed  in  most  literature  under  dimensional  analysis,  the  Bucking- 
ham pi  theorem,  or  the  Rayleigh  theorem.  Dimensional  analysis  can  be 
useful  in  the  development  of  either  distorted  or  undistorted  physical 
models  of  tidal  inlets. 

3.  Purpose  and  Scope  of  Hydraulic  Inlet  Studies. 

Models  of  tidal  inlets  are  constructed  primarily  for  predicting  the 
tidal  currents  and  water  level  fluctuations  in  and  near  inlets.  Predic- 
tions of  the  changes  in  the  flow  regime  associated  with  dredging  or  con- 
struction activities  are  especially  important.  A knowledge  of  currents 
and  tidal  flow  is  needed  for  predicting  sediment  motion  and  the  effects 
of  currents  on  navigation.  Predicting  the  effects  of  currents,  however, 
is  not  discussed  in  this  study. 

Nearshore  currents  and  water  levels  are  often  affected  by  wind  and 
wind- generated  waves  as  well  as  by  astronomical  tides.  The  modification 
of  waves  in  shallow  water  gives  rise  to  radiation  stresses  which  generate 
nearshore  currents,  reduce  the  mean  water  level  in  the  breaker  zone,  and 
increase  the  mean  water  level  on  the  beach.  Wave  motions  also  affect  the 
frictional  effects  on  flows  of  tidal  period.  However,  the  time  and  space 
scales  of  tides  and  surface  waves  are  so  different  that  modeling  the 
details  of  both  at  the  same  time  is  not  feasible.  A study  of  tidal  flows 
requires  dealing  with  no  more  than  the  statistical  effects  of  the  wave 
field.  In  a study  of  wind- generated  waves,  tidal  elevations  and  tidal 
currents  may  be  regarded  as  quasi  steady-state  phenomena.  This  is  also 
true  for  most  tsunamis;  however,  it  may  be  necessary  to  consider  several 
phases  of  the  tide.  In  this  study,  only  the  possible  effects  of  surface 
waves  on  tidal  flows  are  considered. 

Models  of  natural  phenomena  are  used  primarily  for  gathering  essen- 
tial data  at  a convenient  rate  and  at  a lower  cost  than  would  be  possible 
with  the  prototype.  Generally,  this  means  a great  acceleration  of  the 
time  scale,  associated  with  an  increasing  density  of  the  data  gathered 
in  space.  The  model  leads  to  correct  predictions  of  prototype  behavior 
if  the  laws  which  govern  both  model  and  prototype  are  identical.  Models 
which  are  governed  by  simpler  laws  than  the  prototype  can  be  extremely 
useful  if  these  laws  are  nearly  identical  with  those  which  govern  the 
phenomena  of  interest  in  the  prototype,  and  if  phenomena  of  secondary 
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importance  which  lead  to  confusion  at  prototype  scale  are  eliminated  in 
the  model;  e.g.,  the  elimination  of  wind  and  wind-generated  waves  may 
facilitate  the  study  of  tidal  currents.  The  wind  and  wind-generated 
waves  can  always  be  eliminated  from  numerical  models,  and  can  generally 
be  eliminated  from  physical  models. 

Although  the  laws  which  govern  a natural  phenomenon  may  be  fully 
understood  and  expressible  by  a set  of  differential  equations  and  bound- 
ary conditions,  an  exact  solution  of  the  equations  may  be  difficult,  and 
thus  the  construction  of  a physical  model  may  be  desirable.  Mathematical 
models  and  dissimilar  physical  analogs,  such  as  electrical  analogs,  are 
possible  only  to  the  extent  that  the  laws  which  govern  the  phenomenon 
are  known.  When  the  laws  which  govern  a physical  phenomenon  are  under- 
stood well  enough  to  permit  a mathematical  description,  even  physical 
replicas  of  the  actual  phenomena  can  be  regarded  as  analog  computers  for 
solving  a set  of  equations  which  govern  the  model.  Therefore,  the  physi- 
cal as  well  as  the  numerical  models  of  tidal  flow  considered  in  this 
study  are  regarded  as  mathematical  models. 

4.  Preview  of  the  Report. 

Equations  which  govern  the  mean  flow  in  incompressible,  nearly  homo- 
geneous fluid  layers  are  given  in  Section  II,  along  with  the  physical 
interpretation  of  each  term.  The  equations  are  developed  without  ad  hoc 
assumptions  in  Appendix  B,  and  are  divided  into  terms  essential  to  all 
studies  of  tidal  flows — descriptive  terms  of  the  mean  flow  which  are 
essential  to  some  but  not  all  investigations,  and  perturbation  terms 
that  are  usually  neglected.  It  is  often  possible  to  simulate  the  net 
effects  of  the  phenomena  described  by  the  perturbation  terms  with  useful 
accuracy  even  though  the  detailed  structure  of  the  phenomena  is  not  con- 
sidered. The  radiation  stresses  are  normally  neglected  terms  that  account 
for  the  exchange  of  momentum  between  waves  and  currents.  Sager  and  Sea- 
bergh  (in  preparation,  1977)  show  that  the  agreement  between  model  results 
and  prototype  observations  can  be  improved  by  including  the  wave  effect 
in  the  physical  model.  Wave  effects  were  not  included  in  the  original 
plans  for  the  model  studies  and  were  not  considered  in  the  numerical 
models  or  the  calibration  of  the  physical  model. 

General  considerations  for  modeling  tidal  flows  and  their  application 
to  distorted  scale  physical  models  with  particular  reference  to  the  Mason- 
boro  Inlet  model  are  discussed  in  Section  III.  General  features  of  numer- 
ical models  and  their  application  to  two-dimensional  hydrodynamic  models 
such  as  the  Masonboro  Inlet  models  are  discussed  in  Section  IV. 

The  material  in  Sections  III  and  IV  is  not  intended  as  a guide  for  the 
design  and  construction  of  either  physical  or  numerical  models,  but  only 
to  provide  insight  for  understanding  the  Masonboro  Inlet  model  and  other 
similar  models.  Special  features  of  each  model  are  described  in  various 
individual  reports;  however,  the  background  information  needed  for  a 
better  understanding  of  these  reports  has  been  omitted  in  the  descrip- 
tion of  each  model.  This  study  attempts  to  provide  that  information, 
including  a few  references  to  the  material  on  the  individual  models. 
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A surprising  result  of  this  investigation  (see  Sec.  Ill)  is  the 
theoretical  evidence  that  radiation  stresses  can  be  simulated  with  good 
approximation  in  a distorted  scale  physical  model.  This  is  believed  a 
new  theoretical  development.  Sager  and  Seabergh  (in  preparation,  1977) 
report  empirical  evidence  supporting  this  concept  but  do  not  provide  or 
reference  any  theoretical  justification.  They  do  not  claim  the  level  of 
quantitative  accuracy  which  the  development  in  this  study  appears  to 
justify. 

The  effects  of  the  earth's  rotation,  the  surface  wind  field,  and  the 
atmospheric  pressure  gradient  cannot  easily  be  modeled  or  satisfactorily 
simulated  in  the  physical  model  of  Masonboro  Inlet  or  in  most  other 
physical  hydraulic  models.  However,  these  phenomena  are  of  secondary 
importance  in  the  Masonboro  Inlet  study. 

Modeling  or  simulating  the  effects  of  the  surface  wind  field,  the 
atmospheric  pressure  gradient,  and  the  principal  effects  of  the  earth's 
rotation  pose  no  insurmountable  difficulties  for  numerical  models.  Small- 
scale  features  of  the  flow,  with  characteristic  lengths  of  less  than  four- 
or  five-grid  increments,  cannot  be  modeled  with  acceptable  accuracy. 
Features  of  the  flow  with  wavelengths  in  excess  of  10  to  20  space  in- 
crements can  be  modeled  with  satisfactory  accuracy.  Fortunately,  the 
most  important  phenomenon  for  the  Masonboro  model,  the  astronomical  tide 
wave,  has  a length  in  excess  of  20  space  increments. 

In  Section  V,  the  one-dimensional  equations  for  unstratified,  incom- 
pressible flow  in  a channel  of  variable  width  and  depth  are  integrated 
along  the  axis  of  the  channel  to  form  an  equation  relating  flow  in  the 
ocean,  inlet,  and  bay.  The  one-dimensional  flow  equations  are  developed 
from  the  two-dimensional  equations  in  Appendix  C.  For  some  flow  condi- 
tions, it  is  possible  to  combine  all  empirical  coefficients  governing  the 
flow  into  a single  parameter.  This  is  the  essence  of  the  lumped  parameter 
model  introduced  by  Keulegan  (1967).  The  systematic  derivation  given  in 
this  study  facilitates  an  evaluation  of  the  assumptions  involved  in  this 
approach  and  reveals  several  inconsistencies  that  will  not  greatly  de- 
crease the  expected  accuracy  of  the  resulting  calculations  for  Masonboro 
Inlet.  In  applications  to  inlet  problems  in  the  Great  Lakes  where  the 
characteristic  periods  are  much  shorter,  these  inconsistencies  can  lead 
to  significant  errors. 

The  calibration  of  both  physical  and  mathematical  models  is  described 
in  Section  VI.  Calibration,  usually  a trial-and-error  process,  is  neces- 
sary to  compensate  for  the  effect  of  physical  processes  neglected  or  im- 
properly represented  in  the  model,  and  the  lack  of  precise  information 
about  the  prototype. 

The  results  obtained  in  the  various  Masonboro  Inlet  models  are  com- 
pared in  Section  VII;  a summary  with  conclusions  and  recommendations  is 
presented  in  Section  VIII. 
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II.  A MATHEMATICAL  DESCRIPTION  OF  TIDAL  FLOWS 


1.  General . 

Tidal  currents  are  quasi-horizontal  and  almost  periodic  with  promi- 
nent components  of  approximately  12  and  24  hours.  Tidal  components  with 
periods  as  short  as  3 hours  can  be  identified  at  some  locations  by  a 
harmonic  analysis  of  the  tide  records.  Although  the  tides  are  affected 
by  the  earth's  rotation,  this  effect  is  insignificant  in  many  estuaries. 
Friction  is  a prominent  factor  in  shallow  areas;  however,  the  effect  of 
friction  is  negligible  in  deep  water.  Most  tidal  flows  are  adequately 
described  by  the  linearized  hydrodynamic  equations.  However,  some  signs 
of  nonlinear  behavior  of  the  tides  in  estuaries,  such  as  overtides  and 
tidal  bores,  are  not  described  adequately  by  the  linearized  equations. 

Meteorological  disturbances  produce  transient  currents  and  water 
level  disturbances  with  time  scales  which  overlap  those  of  the  astro- 
nomical tides.  A precise  separation  of  the  meteorological  and  astro- 
nomical effects  on  sea  level  and  nearshore  currents  is  never  easy  and  is 
not  always  possible.  The  overlap  occurs  in  the  analysis  of  both  empirical 
data  and  the  hydrodynamic  equations.  Thus,  any  thorough  study  of  the  flow 
in  tidal  inlets  and  estuaries  must  consider  both  meteorological  and  astro- 
nomical effects.  The  effects  of  the  variation  of  gravitational  attraction 
of  the  sun  and  moon  over  a typical  estuary  are  so  small  that  the  tides  in 
most  estuaries  are  adequately  represented  by  prescribing  appropriate  bound- 
ary conditions  in  the  open  sea. 

Many  distinct  hydrodynamic  phenomena  are  associated  with  the  tides, 
but  not  all  are  important  at  any  one  location.  Some  phenomena  are  promi- 
nent in  parts  of  a problem  area  or  for  short-time  periods,  and  insignifi- 
cant in  other  parts  of  the  area  and  at  other  times.  Thus , any  set  of 
equations  which  provides  a satisfactory  description  of  a wide  range  of 
tidal  problems  must  be  fairly  complete  and  must  contain  terms  which  can 
be  neglected  in  many  specific  applications. 

The  most  reliable  method  of  formulating  a dependable  set  of  equations 
which  satisfies  the  above  requirements  is  to  begin  with  a complete  set  of 
hydrodynamic  equations  which  describe  most  fluid  flows  and  systematically 
modify  them  to  permit  a ready  separation  of  the  phenomena  of  primary 
interest  from  phenomena  of  generally  secondary  importance. 

A derivation  of  the  long  wave  equations  (carried  out  as  discussed 
above)  is  given  in  Appendix  B.  An  identical  set  of  equations  for  the 
terms  of  primary  importance  has  been  derived  directly  by  many  authors  by 
simply  neglecting  all  phenomena  considered  of  se con ’ary  importance.  The 
most  important  terms  have  been  obtained  by  establishing  a balance  of 
forces  on  a fluid  element  by  considering  only  the  obviously  important 
terms.  A rigorous  derivation  is  used  in  this  study  to  provide  at  least 
subjective  estimates  of  the  errors  inherent  in  the  approximating  equa- 
tions and  to  show  that  the  errors  committed  in  numerical  and  physical 
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modeling  differ.  Sometimes  one  type  of  error  is  more  acceptable  than 
another  in  solving  a specific  problem. 

Sound  waves  are  eliminated  from  the  fundamental  equation  by  assuming 
at  the  outset  that  the  water  is  incompressible.  The  detailed  structure 
of  wind-generated  waves  and  turbulence  is  eliminated  from  the  equations 
without  eliminating  the  interaction  effects  by  averaging  each  important 
variable  over  a finite  volume  and  time  interval. 


The  practical  necessity  of  eliminating  some  detail  of  the  flow  may 
be  recognized  when  it  is  realized  that  all  physical  measurements  of 
velocity  are  averages  over  some  finite  volume  and  time  period,  if  only 
because  of  the  physical  size  and  finite  response  time  of  the  measurement 
instruments.  Displacements  of  the  fluid  smaller  than  the  effective  size 
of  the  instruments,  and  motions  with  time  scales  smaller  than  the  response 
times  of  the  instruments  cannot  be  accurately  measured.  The  effective 
averaging  volume  and  the  time  interval  referred  to  prototype  scale  are 
generally  larger  for  physical  models  than  for  the  prototype,  and  larger 
for  numerical  models  than  for  physical  models.  Recognition  of  this  dif- 
ference is  essential  to  the  proper  interpretation  of  the  results  of  model 
studies. 

Details  of  the  vertical  structure  are  eliminated  by  integrating  the 
flow  from  the  seabed  to  the  free  surface.  Thus  the  horizontal  velocities 
are  replaced  by  total  fluid  fluxes,  or  rate  of  transport,  and  the  verti- 
cal velocity  is  represented  only  by  the  stretching  or  compression  of 
fluid  layers.  This  mode  of  derivation  leads  to  a set  of  equations  which 
approximately  describe  both  tidal  flows  and  storm  surges,  and  a set  of 
remainder  terms  which  can  be  used  to  estimate  the  error  committed  in 
using  the  approximate  equations. 


2.  Motions  of  Tidal  Scale. 

The  equations  for  the  conservation  of  momentum  in  quasi-horizontal 
flow  (App.  B,  eqs.  B39  and  B40)  may  be  expressed  in  the  form: 


(a)  (b) 


(c) 


(d) 


(e) 
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The  two-dimensional  equation  of  continuity  (App.  B,  eq.  B23)  is: 


(h)  (i)  (j) 

k<fi  - D>  * I?  * if  ■ R * E ■ (3) 


where 

t 


= time 


x,  y , z 

E 


U 


V 

R 


P 


= principal  axes  in  a right-hand  coordinate  system; 
z is  positive-upward 


h 

= J udz,  mean  fluid  transport  parallel  to  the  x-axis 
Z 

h 

= J vdz,  mean  fluid  transport  parallel  to  the  y-axis 
Z 

= precipitation  on  the  water  surface  and  bottom 
seepage  into  the  basin  in  excess  of  evaporation, 
and  bottom  seepage  out  of  the  basin 

h 

= F / pdz 

Z 


= density 


g = acceleration  of  gravity 

u,  v = mean  horizontal  velocity  components  parallel  to  the 

x-  and  y-axes,  where  the  mean  is  taken  over  a small 
time  and  space  interval 

h = h(x,y,t),  mean  elevation  of  the  top  of  a fluid  layer 

above  datum 

D = D(x,y,t)  = h - Z,  total  depth 

Z = Z(x,y,t),  mean  elevation  of  the  bottom  of  a fluid 

layer  above  datum 

f = sin  , Coriolis  parameter 
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= angular  speed  of  the  earth's  rotation, 
7.292  x 10" 5 radians  per  second 

= latitude 


? = p(x,y,z,t),  pressure  within  the  fluid 

p^  = pressure  at  the  top  of  a fluid  layer 

(t zx^h’  ^zy^h  = surface  stresses  parallel  to  the  x-axis  at  the  top 
of  a fluid  layer 

(Tzy^Z  = surface  stresses  parallel  to  the  x-  and  y-axes  at 
the  bottom  of  a fluid  layer 

u = u(x,y,z,t)  and  v = v(x,y,z,t),  horizontal  components 

of  fluid  velocities  parallel  to  the  x-  and  y-axes 

S;cw>  Syw>  = interaction  between  the  mean  flow  and  the  perturb a- 

Sxt»  tions  due  to  waves  or  turbulence 

$xo  > Syo > 

S xv  > $yv>  E 

The  contribution  of  rainfall,  evaporation,  and  seepage  to  the  momen- 
tum equations  has  been  omitted  because  these  are  expected  to  be  smaller 
than  the  uncertainty  in  evaluating  the  stress  terms  in  most  applications; 
e.g.,  the  inlet  problems  considered  in  this  study. 


Essential  terms  for  the  discussion  of  all  two-dimensional  tidal  flows 
in  shallow  water  are  shown  on  the  left  of  the  equal  sign  in  equations  (1), 
(2),  and  (3).  Important  terms  for  some  but  not  all  flows  of  this  class 
are  shown  on  the  right. 


3.  Interpretation  of  the  Principal  Terms. 

The  lower-case  letters  in  parenthesis,  (a)  to  (j)  , over  the  terms  in 
the  equations  are  used  to  identify  the  terms  in  the  following  discussion 
of  their  physical  significance  with  some  common  approximations. 

(a)  3U/3t  and  3V/3t.  Identifies  the  change  of  the  integrated  flow 
per  unit  width  with  time.  This  change  may  result  from  the  temporal  accel- 
eration of  the  mean  flow  or  a change  in  thickness  of  the  fluid  layer.  At 
least  one  of  the  two  terms  is  important  in  every  tidal  or  storm  surge 
problem. 

(b)  gD3h/3x  and  gD3h/3y.  Describes  the  slope  of  the  water  surface. 
The  pressure  gradients  which  result  from  the  sloping  water  surface  are 
the  principal  driving  forces  for  all  purely  tidal  flows  in  estuaries. 

They  are  also  prominent  in  all  storm  surge  problems. 

(c)  (Tzx)z  an< 1 (T3j/)z*  Usually  interpreted  as  the  stress  of  the 
fluid  layer  against  the*DOttom  boundary,  but  in  many  applications. 
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internal  stresses  are  implicitly  included.  In  most  problems,  the  bottom 
boundary  is  the  solid  earth  and  is  assumed  impermeable.  These  terms  are 
usually  treated  as  the  bottom  friction  terms.  In  studying  stratified 
flow,  it  may  be  necessary  to  recognize  two  distinct  fluid  layers.  When 
this  is  true,  (x3X)z  and  (Tzi/^Z  f°r  the  uPPer  layer  may  be  the  stress 
of  the  saline  water  in  the  lower  layer  on  the  freshwater  in  the  upper 
layer. 

The  stress  within  a fluid  has  the  dimension  of  the  product  of  density 
and  the  square  of  the  speed.  In  general,  the  magnitude  of  the  stress 
terms,  (t 3x)z  and  (t zy)z,  increases  with  the  current  speed.  Thus,  these 
terms  are  naturally  represented  by  the  expression: 

(”)  xss  = Kq2  , (4) 

where  K is  a dimensionless  coefficient  which  may  be  a function  of  speed 
and  other  parameters  governing  the  flow,  q is  the  velocity  parallel  to 
the  boundary,  and  the  s-axis  is  chosen  parallel  to  the  flow.  The  bottom 
stress  vector  is  generally  assumed  parallel  to  the  velocity,  q,  which 
generally  increases  with  distance  from  the  nearest  boundary.  In  speci- 
fying a horizontal  surface  stress,  K must  be  a function  of  the  distance 
between  the  boundary  and  the  point  at  which  q is  specified.  Most  hydro- 
dynamics and  hydraulics  textbooks  show  that  K increases  as  the  roughness 
of  the  surface  increases. 

In  a study  of  quasi-horizontal  flows,  the  velocity  profile  as  a func- 
tion of  distance  above  the  bottom  boundary  of  the  fluid  is  generally  not 
available.  Therefore,  the  best  available  approximation  is  usually  of  the 
form: 

(j)(TZS)z  = KlQlgy  , (5) 

where  K decreases  as  the  depth  increases  and  increases  as  the  bottom 
roughness  increases,  and  Q = (U2  + V2)1/2  and  has  the  same  direction  as 
the  mean  flow.  The  stresses,  CTsar) Z and  ixzy^Z>  are  the  components  of 
(xzs)Z  parallel  to  the  x-  and  y-axes . 

The  coefficient,  K,  may  be  approximated  by  the  expression: 

K = K2D"a2  , (6) 

where  K2  is  a parameter  used  to  characterize  the  bottom  roughness  and 
d2  is  an  empirical  parameter.  Manning's  investigation  of  expressions 
similar  to  equation  (6)  was  summarized  by  Chow  (1959,  p.  100).  Initially, 
a 2 was  regarded  as  a constant  less  than  one.  Ultimately,  it  was  decided 
that 
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The  value  1/3  has  been  adopted  in  most  studies;  however,  a few  continue 
to  use  other  values  of  the  exponent.  When  the  metric  system  is  used 
with  Hanning's  formulation,  the  bottom-stress  coefficient  for  use  in 
equation  (5)  is  given  by: 


K = 


sn! D 1 / 3 

(m/s3)2/3 


(7a) 


where  n is  a dimensionless  constant  which  characterizes  the  bottom 
roughness.  When  the  U.S.  Customary  system  is  used. 


SH 1 D 1 / 3 

(3.2808  ft/s3)2/3 


(7b) 


The  bottom  roughness  depends  on  the  bottom  bed  forms.  The  bed  forms 
in  sedimentary  materials  vary  with  the  current  speed  and  the  composition 
of  the  bed.  In  general,  the  effective  roughness  increases  with  the 
current  speed  to  some  small  finite  speed  and  decreases  for  higher  speeds. 
Ludwick  (1974)  reported  620  measurements  of  the  bottom-stress  coefficients 
near  the  entrance  to  Chesapeake  Bay,  and  found  that  K varied  widely  with 
the  development  of  bed  forms.  In  Ludwick 's  data,  the  characteristic  value 
of  K for  ebb  flows  was  about  three  times  the  characteristic  value  for 
floodflows . 


Behrens,  Watson,  and  Mason  (1977)  measured  n in  a Texas  tidal 
inlet  and  also  found  that  the  bottom  roughness  varied  with  the  develop- 
ment of  bed  forms,  but  did  not  find  a significant  difference  between  ebb 
and  floodflows.  Odd  (1971)  reported  that  the  friction  factor,  K,  also 
varied  with  the  phase  of  the  tide  even  with  no  change  in  the  bottom  rough- 
ness, but  he  does  not  provide  any  quantitative  information  or  references 
to  other  sources  where  quantitative  information  may  be  found. 

From  the  above  discussion  it  seems  that  the  true  bottom  friction 
depends  primarily  on  the  velocity  near  the  bottom  and  the  small-scale 
features  of  the  flow  that  are  not  considered  in  detail.  Several  phe- 
nomena occur  in  nature  to  produce  bottom  stresses  which  cannot  be 
described  by  the  quadratic  stress  law;  e.g.,  a finite  time  is  required 
for  any  contribution  of  the  surface  wind  to  the  current  to  be  experienced 
at  the  bottom.  The  time  required  increases  with  increasing  water  depth 
and  the  direction  of  the  bottom  velocity  may  differ  from  that  of  the 
mean  velocity.  If  the  wind  over  a small  basin  continues  long  enough  for 
an  equilibrium  slope  to  develop,  the  mean  velocity  may  vanish  and  the 
bottom  velocity  may  be  directed  opposite  to  the  free-surface  velocity. 
Other  phenomena  which  lead  to  differences  between  the  mean  velocity  and 
the  bottom  velocity  are  discussed  later. 


Thus,  when  equations  (4)  or  (5)  are  used  to  evaluate  the  bottom 
stress,  K should  increase  with  increasing  roughness  of  the  bottom  and 
decrease  with  increasing  depth.  K may  also  increase  with  increasing 
windspeed  or  wave  height. 
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Bottom  friction  dissipates  energy  and  thus  contributes  toward  a 
reduction  in  the  amplitude  of  the  tide  wave  and  a decrease  in  the  phase 
velocity  of  gravity  waves. 

(d)  anc*  ( Tzy^h'  Usually  interpreted  as  the  stress  in  the 

x-  and  y-directions  of  an  upper  fluid  layer  on  a lower  fluid  layer.  When 
two  distinct  fluid  layers  are  identified,  and  [tzy)h  for  the  lower 

layer  may  be  recognized  as  the  stress  exerted  on  the  lower  layer  by  the 
upper  layer. 

At  a free  water  surface,  these  terms  describe  the  friction  between 
the  wind  and  the  water.  The  direct  effect  of  the  wind  stress  is  to  gen- 
erate a surface  current  in  the  direction  of  the  wind.  If  the  movement 
of  this  current  is  impeded  by  a barrier,  the  current  will  be  partially 
replaced  by  a slope  of  the  water  surface,  upward  in  the  direction  of  the 
airflow.  The  wind  stress  is  an  essential  term  in  the  equations  for  storm 
surge,  and  is  generally  neglected  in  investigations  of  tidal  flows  unless 
believed  an  essential  element  of  the  problem.  However,  wind  usually  has 
some  effect  on  the  water  level  and  currents  in  all  tidal  basins.  Neglect- 
ing the  wind  when  it  is  a significant,  though  perhaps  secondary  factor, 
may  lead  to  errors  in  calibrating  and  verifying  a model  for  purely 
gravity  flows . 

(e)  9/9x  U2/D  + 9/ 9y  UV/D  and  9/9x  UV/D  + 9/9y  V2/D.  May  be 
interpreted  as  integrated  measures  of  the  advection  of  momentum  or  inte- 
grated measures  of  the  inertia  of  the  water.  These  terms  may  be  impor- 
tant only  in  regions  with  strong  velocity  gradients,  such  as  the  entrances 
and  exits  of  inlets,  junctions  of  canals,  flood  plains,  or  other  regions 
where  the  flow  is  compelled  to  undergo  large  changes  in  direction. 

(f)  fU  and  fV.  The  Coriolis  terms  which  result  from  the  neces-  -v. 

sity  of  using  a coordinate  frame,  fixed  to  the  rotating  earth;  Newton's 

laws  of  motion  apply  to  a coordinate  system  fixed  with  respect  to  the 
stars.  The  Coriolis  effect  is  discussed  here  in  detail  and  documented, 
because  the  conditions  under  which  the  Coriolis  effects  are  important 
are  frequently  misinterpreted  in  discussions  of  tidal  hydraulics. 

The  earth's  rotation  can  affect  the  water  motion  in  several  distinct 
ways,  none  of  which  are  considered  significant  in  typical  inlets  with 
widths  of  about  1 mile  or  less  and  maximum  currents  of  about  2 meters 
(6.6  feet)  per  second  or  less.  The  Coriolis  effect  is  measurable  in 
some  of  the  wider  and  deeper  inlets,  or  when  current  speeds  are  excep- 
tionally high.  Coriolis  acceleration  is  often  a prominent  feature  of 
the  tidal  flow  in  the  open  sea  adjacent  to  the  inlet  and  it  may  have  a 
significant  effect  in  the  bay  or  lagoon  on  the  landward  side  of  the 
inlet.  Flow  resulting  from  the  Coriolis  acceleration  is  often  called 
the  inertial  flow. 

The  earth's  rotation  tends  to  produce  an  apparent  rotation  of  the 
fluid  toward  the  right  in  the  Northern  Hemisphere  and  toward  the  left 
in  the  Southern  Hemisphere.  The  period  of  this  rotation,  called  the 
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inertial  period,  is  equal  to  1/f.  The  circle  traversed  by  a water  par- 
ticle is  called  the  inertial  circle,  its  radius,  Rg,  is  given  by: 
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_ v 

" 7 ’ 


(8) 


where  v is  the  current  speed. 

When  the  flow  is  generated  by  pressure  gradients  in  a homogeneous 
fluid,  as  for  tidal  flows  near  the  coast,  the  direction  and  magnitude 
of  the  inertial  flow  are  nearly  independent  of  depth.  Verber  (1966) 
reported  observations  of  the  inertial  flow  in  Lake  Michigan  with  currents 
as  low  as  6 centimeters  (0.2  foot)  per  second.  The  radius  of  the  corres- 
ponding inertial  circle  is  about  586  meters  (1,923  feet). 

When  the  flow  is  generated  by  frictional  stresses  at.  the  fluid  bound- 
aries, as  for  wind  stress  at  the  free  surface  and  bottom  friction,  a 
finite  time  is  required  for  the  rotational  effect  to  penetrate  to  the 
interior  of  the  fluid.  As  a result,  there  is  an  additional  rotation  to 
the  right  with  distance  from  the  fluid  boundary.  The  frictional  force 
at  the  bottom  is  directly  opposed  to  the  bottom  current  which  would  exist 
in  the  absence  of  friction.  These  influences  cause  the  vertically  inte- 
grated flow  vector,  resulting  from  frictional  forces,  to  be  directed  to 
the  right  of  both  the  surface  and  bottom-stress  vectors.  The  fluid 
layers  affected  by  this  rotating  current  are  called  the  Ehvan  boundary 
layers.  The  thickness  of  the  layers  increases  with  the  intensity  of 
turbulence  in  the  fluid.  Thus,  for  turbulent  flow  in  shallow  water,  the 
two  boundary  layers  may  overlap  and  little  rotation  of  the  currents  is 
observed.  The  combined  effects  of  the  earth's  rotation  and  fluid  fric- 
tion on  wind-generated  and  tidal  currents  are  discussed  by  Defant  (1961; 
Vol.  1,  Chs.  10  and  13;  Vol.  2,  Ch.  1)  and  Neuman  and  Pierson  (1966; 

Ch.  8).  Defant  also  discussed  the  inertial  circle.  Webster  (1969)  sum- 
marized most  observations  of  the  inertial  circle  published  before  1968. 
Platzman  (1963),  Jelesnianski  (1967,  1970),  and  Forristall  (1974)  pre- 
sented schemes  for  including  the  effect  of  the  earth's  rotation  on  the 
stress  vectors  in  storm  surge  computations.  If  this  rotation  is  in- 
hibited by  lateral  boundaries,  it  is  replaced  by  an  upward  slope  to  the 
right  of  the  current  in  the  Northern  Hemisphere  and  to  the  left  in  the 
Southern  Hemisphere. 

The  magnitude  of  this  slope  can  be  computed  from  equation  (2)  by 
considering  only  the  slope  and  Coriolis  terms  to  obtain: 
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where  the  y-axis  is  assumed  to  coincide  with  the  axis  of  the  inlet  and 
v is  the  mean  velocity  of  the  flow,  parallel  to  the  y-axis  averaged  in 
the  vertical.  Thus,  the  slope  across  the  inlet  is  given  by: 


Ah 


(10) 


In  applying  the  above  considerations  to  Masonboro  Inlet,  field  meas- 
urements show  that  the  inlet  is  well  mixed.  This  implies  that  the  top 
and  bottom  Ekman  layers  overlap  so  that  little  rotation  of  the  current 
vector  with  depth  is  expected  to  result  from  the  rotation  of  the  earth. 

The  maximum  current  speeds  recorded  in  Masonboro  Inlet  are  about  1 meter 
(3.3  feet)  per  second.  The  Coriolis  parameter,  f,  for  latitude  34°  is 
8.155  x 10  5 per  second,  and  the  width  of  the  inlet  is  457  meters  (1,500 
feet).  Since  the  radius  of  the  inertial  circle  is  about  12.5  kilometers 
(7.8  miles)  or  R = v/f,  the  inertial  circle  cannot  develop.  The  net 
change  in  elevation  across  the  narrow  inlet,  as  computed  by  equation  (10), 
is  less  than  4 millimeters  (0.16  inch);  however,  this  effect  is  insignifi- 
cant. Reid,  Vastano,  and  Reid  (in  preparation,  1977)  show  that  bottom 
friction  in  shallow  areas  may  overshadow  the  Coriolis  acceleration  during 
storms  when  current  velocities  are  high,  even  though  Coriolis  accelera- 
tion and  bottom  friction  are  of  the  same  magnitude  during  low  flow  con- 
ditions . 


(g)  (l/p)(D)  3p^/3x  and  (1/p) (D)  3p^/3y.  Describes  the  com- 
ponent of  the  atmospheric  pressure  gradient  which  affects  the  water 
motion.  The  atmospheric  pressure  gradient  is  generally  neglected  in 
studies  of  water  motions  unless  believed  an  essential  factor  in  the 
phenomena  being  studied.  The  atmospheric  pressure  gradient  terms  are 
generally  of  trivial  importance  in  studies  of  tidal  flows,  except  for 
hurricane-generated  storm  surges. 

(h)  3 (h  - Z)/3t.  A measure  of  the  time  rate  of  expansion  or 
contraction  of  a fluid  layer.  A well-mixed  estuary,  such  as  Masonboro 
Inlet,  can  be  adequately  represented  as  a single  fluid  layer.  In  this 
case,  h is  the  elevation  of  the  free  surface  and  Z is  the  elevation 
of  the  seabed.  Later  in  this  study,  Z is  assumed  independent  of  time. 
With  this  assumption,  the  term  3h/3t  in  equation  (3)  measures  the  rate 
at  which  the  water  level  is  rising  or  falling. 

(i)  3U/3x  + 3V/3y.  A measure  of  the  rate  at  which  the  water  is 
converging  or  diverging  horizontally  at  any  point  (x,  y) . The  complete 
continuity  equation  (3)  states  that  the  water  level  must  rise  in  an  area 
of  horizontal  convergence  and  fall  in  an  area  of  horizontal  divergence. 
At  least  one  of  these  terms  is  essential  to  all  tide  and  storm  surge 
models . 


(j)  R0.  A measure  of  the  precipitation  on  the  water  surface  and 

bottom  seepage  into  the  basin  in  excess  of  the  evaporation  and  bottom 
seepage  out  of  the  basin.  In  general,  R is  not  important  in  tidal 
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hydraulics  investigations.  The  most  prominent  exception  is  a nearly 
landlocked  estuary  or  lagoon  during  severe  rainstorms. 

4.  The  Perturbation  Terms. 


The  principal  terms  in  equations  (1),  (2),  and  (3)  were  derived  by 
the  vertical  integration  and  time  averaging  of  the  Navier-Stokes  equa- 
tions. Deviations  between  the  actual  flows  and  the  flows  defined  by 
vertical  integration  and  time  averaging  may  be  defined  as  perturbations. 
Nonlinear  terms  in  the  perturbations  represent  the  interaction  between 
the  small-scale  features  of  the  flow,  such  as  turbulence,  surface  waves, 
and  the  vertical  shear  of  the  current.  These  phenomenon  are  discussed 
below.  The  turbulent  stresses  (App.  B,  eqs.  B41  and  B42)  are: 
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The  Reynolds  stress  terms,  (u1)2,  (v1)2,  and  u'v1 , are  generally  much 
larger  than  the  terms  involving  molecular  viscosity  except  in  the  thin 
viscous  sublayer,  never  more  than  a few  millimeters  in  thickness.  The 
viscous  terms  can  generally  be  neglected  when  discussing  prototype  flows. 
The  simple  form  of  the  expression  describing  viscous  shears  can  be  re- 
covered by  introducing  mixing  coefficients  defined  by  the  relations : 
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The  mixing  coefficients  will  be  needed  later  in  the  discussion  of  numeri- 
cal models. 


The  mixing  coefficients,  sometimes  called  macroviscosity  or  "eddy" 
coefficients,  depend  on  the  structure  of  the  flou  field.  In  contrast, 
the  viscosity  coefficient  depends  on  the  -physical  properties  of  the 
fluid.  The  mixing  coefficients  cannot  be  evaluated  from  theoretical 
considerations,  and  are  usually  determined  empirically  from  measured 
stresses  and  velocity  gradients,  or  selected  subjectively  to  obtain  an 
acceptable  solution  to  a problem.  Fortunately,  the  solution  of  many 
tidal  hydraulics  problems  involving  the  mixing  coefficients  is  not 
highly  sensitive  to  the  values  selected. 
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5.  The  Radiation  Stresses. 


The  earliest  investigations  of  the  effects  of  the  Reynolds  stresses 
due  to  waves  on  the  mean  flow  were  reported  by  Longuet-Higgins  and 
Stewart  (1962,  1963,  1964).  The  contribution  of  the  radiation  stresses 
to  the  mean  transport  equations  as  developed  in  Appendix  B (eqs.  B43  and 
B44)  is  given  below  in  the  notation  of  Longuet-Higgins  and  Stewart. 
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where  u,  v,  and  w are  the  components  of  motion  due  to  waves,  and  h 
is  the  instantaneous  elevation  of  the  free  surface. 

The  interpretation  of  these  stresses  may  be  clarified  by  considering 
a train  of  gravity  waves  traveling  in  the  s-direction,  at  an  angle  6 to 
the  x-axis.  q is  the  horizontal  particle  velocity  in  the  s-direction, 
due  to  the  gravity  waves.  For  these  conditions. 
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The  theory  of  gravity  waves  shows  that  q2,  w2,  and  fT2  are  propor- 
tional to  A2  where  A is  the  amplitude  of  the  wave.  Equations  (1), 

(2) , (19) , and  (20)  show  that  a decrease  in  wave  amplitude  in  the  direc- 
tion of  wave  propagation  is  associated  with  an  increase  in  the  momentum 
of  the  mean  flow  or  an  upward  slope  of  the  water  surface  in  the  direction 
of  wave  propagation  and  vice  versa. 

The  computed  magnitude  of  the  radiation  stresses  depends  on  the 
assumed  wave  shape.  Longuet-Higgins  and  Stewart  (1962,  1963,  1964) 
used  the  theory  of  infinitesimal  waves  to  provide  estimates  for  the 
radiation  stress  terms.  James  (1974a,  1974b)  used  cnoidal  and  hyper- 
bolic wave  theory  to  provide  similar  estimates. 

The  component  of  radiation  stress  parallel  to  the  beach  is  responsi- 
ble for  the  generation  of  the  longshore  currents,  important  in  sediment 
transport.  The  component  of  radiation  stress  normal  to  the  beach  is 
responsible  for  wave  setup;  i.e.,  an  upward  tilt  of  the  water  surface 
from  the  breaker  zone  to  the  beach . 

6.  Other  Secondary  Terms. 

A second-order  contribution  to  the  continuity  equation  resulting  from 
wave  action  is  given  by  equation  (B25)  in  Appendix  B as : 

E = - [|^h  ♦ f^h]  , (21) 

where  the  subscript  h indicates  that  u and  v are  to  be  evaluated 
at  z = h. 

This  term  represents  the  fluid  transport  between  the  trough  and  crest 
of  the  waves,  and  is  analogous  to  the  radiation  stresses  discussed  above. 
The  term  can  be  significant  only  when  the  wave  height  is  of  the  same 
order  of  magnitude  as  the  water  depth  and  is  changing  rapidly  in  space. 
Thus,  the  term  may  make  a significant  contribution  to  the  conservation 
of  fluid  volume,  but  only  in  and  near  the  surf  zone  and  over  shallow 
shoals . 

A small  error  is  made  in  using  the  product  of  the  integrals  of  the 
horizontal  velocity  components  in  place  of  the  integral  of  the  products. 
These  errors,  represented  by  and  S yV,  have  the  form: 
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where 


U V 

u"  = u - -jj  , v"  = v - p . (24) 

The  terms  u"  and  v"  are  the  deviations  of  the  true  mean  velocity 
components,  as  functions  of  z from  the  vertically  averaged  velocity. 
Derivatives  of  these  integrals  will  be  small  except  in  regions  of  rela- 
tively steep  changes  in  bottom  topography.  These  terms  contribute  toward 
a downward  slope  in  the  water  surface  or  a deceleration  of  the  flow  on  an 
upward  slope  and  a change  of  opposite  sign  on  a downward  slope.  Where 
the  depth  change  is  due  to  a submerged  bar,  the  change  in  flow  character- 
istics is  reversible.  Neglecting  these  terms  in  calculations  will  lead 
to  an  error  in  computed  results  near  the  bar,  but  may  not  have  a detect- 
able effect'  a short  distance  away. 


Perturbation  terms  arising  from  products  of  the  mean  flow  and  pertur- 
bation variables,  or  from  time  derivatives  of  the  perturbation  variables, 
are  given  by: 


where  all  velocity  terms  are  evaluated  at  z = h. 

If  the  bottom  contours  do  not  change  noticeably  within  a single  wave 
period. 


ws  = uZa7  + vz^r  • C27) 

Thus,  for  a single  wave  train,  the  final  term  in  equations  (25)  and  (26) 
has  the  form: 
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where  q is  the  horizontal  particle  velocity  near  the  bottom  due  to 

Lj 

wave  motions,  (8Z/3s)  is  the  slope  of  the  bottom  in  the  direction  of 
maximum  slope,  and  0 is  the  angle  between  the  direction  of  wave  propa- 
gation and  the  direction  of  maximum  slope.  Mean  slopes  large  enough  to 
produce  measurable  effects  cannot  exist  over  any  appreciable  distance. 

The  remaining  terms  in  equations  (25)  and  (26)  describe  interactions 
between  surface  waves  and  other  features  of  the  flow.  These  terms  have 
not  been  investigated  in  detail,  but  they  can  make  a significant  contri- 
bution to  the  momentum  balance  only  in  a region  where  the  wave  height  is 
a significant  fraction  of  the  mean  depth.  The  contribution  is  greatest 
when  the  wave  height  is  changing  rapidly  with  time  or  distance.  The 
conditions  which  favor  the  importance  of  these  terms  are  most  likely 
to  exist  in  the  surf  zone. 

7.  Boundary  Conditions. 

In  developing  the  two-dimensional  equations  of  motion,  it  has  been 
assumed  that  no  fluid  crosses  upper  and  lower  boundaries.  The  possi- 
bility of  rainfall  or  evaporation  was  treated  by  modifying  the  continuity 
equation.  The  conditions  of  the  flow  at  lateral  boundaries  must  be  speci- 
fied before  the  two-dimensional  equations  can  be  solved.  Three  types  of 
boundary  conditions  may  occur  when: 

(a)  The  flow  across  a fixed  boundary  is  specified  as  a 
function  of  time  along  the  boundary; 

(b)  no  flow  occurs  across  parts  of  the  boundary,  but  the 
location  of  the  boundary  varies  in  time,  and  determining  the 
location  of  the  boundary  is  part  of  the  problem;  and 

(c)  the  elevation  of  the  water  surface  is  specified  as  a 
function  of  time. 

The  application  of  boundary  conditions  is  an  important  part  of  model 
design  and  is  discussed  in  the  following  sections. 

III.  MODELING  TIDAL  FLOWS:  PHYSICAL  MODELS 


1.  General. 


Many  characteristics  of  tidal  waterways  can  be  simulated  with  useful 
accuracy  without  an  exact  reproduction  of  all  the  processes  described 
previously.  An  advantage  of  modeling  is  the  ability  to  eliminate  some 
of  the  secondary  factors  from  a study  of  phenomena  of  primary  importance. 

Physical  models  of  tidal  flows  are  generally  constructed  as  small- 
scale,  geometrically  similar  replicas  of  the  prototype  with  the  excep- 
tions that  the  saltwater  of  the  prototype  (if  present)  is  usually 
replaced  by  freshwater  and  the  scale  ratio  used  for  the  vertical 
dimension  may  differ  from  that  used  for  horizontal  dimensions. 
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Rouse  and  Ince  (1957)  reported  that  Reynolds  in  1887  showed  that  for 
model  studies  which  do  not  involve  intermediate  or  deepwater  waves,  the 
vertical  and  horizontal  scales  for  tidal  flows  could  be  selected  inde- 
pendently but  that  the  time  scale  must  be  a function  of  the  vertical  and 
horizontal  space  scales.  When  the  horizontal  and  vertical  scales  are 
selected  independently,  the  model  is  considered  distorted.  The  undis- 
torted model  may  be  regarded  as  a special  case  of  the  distorted  model  in 
which  the  ratio  of  the  vertical  and  horizontal  scales  is  unity. 

Some  terms  in  the  equations  of  motion  (1),  (2),  and  (3),  e.g.,  gravity, 
cannot  be  scaled  in  the  laboratory.  Others,  such  as  the  rotation  of  the 
earth,  the  atmospheric  pressure,  and  wind-shear  forces  on  the  water,  can 
be  modeled  but  only  with  great  difficulty  and  expense. 

These  facts  must  be  carefully  considered  in  designing  the  model. 
Although  it  is  impossible  to  properly  scale  a wide  range  of  physical 
processes  at  one  time,  it  is  often  possible  to  scale  the  most  important 
processes  for  a given  phenomenon  and  to  simulate  the  effects  of  many 
secondary  processes  on  the  flows  of  primary  importance. 

Since  these  simulated  motions  do  not  model  the  prototype  flow,  the 
motions  cannot  be  controlled  correctly  from  theoretical  considerations 
alone  and  some  adjustment  of  the  model  is  usually  required  to  obtain 
optimum  results.  This  adjustment  process  is  included  in  model  calibra- 
tion. 

Physical  models  have  the  distinct  advantage  of  permitting  a modeling 
of  the  boundaries  of  flow  as  accurately  as  desired  without  significantly 
increasing  the  operating  cost  of  the  model.  Physical  models  are  also 
useful  for  visually  demonstrating  many  phenomena  of  potential  importance. 

Both  physical  and  numerical  models  are  useful  in  understanding 
natural  phenomena.  An  account  of  the  advantages  and  disadvantages  of 
each  model  is  as  important  an  aspect  of  this  evaluation  as  a comparison 
of  the  results  achieved  in  several  experiments.  An  examination  of  the 
principal  features  of  each  type  of  model  is  necessary  to  proceed  further 
with  the  comparison.  Physical  models  are  discussed  in  this  section; 
numerical  models  are  discussed  in  Section  IV. 

2.  Physical  Models  for  Quasi-Horizontal  Flow. 

The  physical  model  used  for  the  study  of  Masonboro  Inlet,  like  most 
tidal  hydraulic  models,  is  a distorted  scale  model. 

The  relation  between  the  scaling  factors  for  time,  space,  and  bottom 
stress  may  be  derived  systematically  by  developing  the  governing  equations 
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in  dimensionless  form,  and  then  comparing  the  dimensional  realizations 
of  these  equations  for  prototype  and  model  scales.  Let 


= «px* 

^ 

= apy* 

= BpZ* 

zm  = V* 

= Ypt* 

~ Y^t*  > 

(29) 

where 

P 


= prototype 
= model 

= dimensionless  variable 


a , am  = characteristic  horizontal  dimensions  of  the  prototype  and 
model 

Bp,  Bm  = characteristic  vertical  dimensions 

Yp,  ym  = characteristic  time  dimensions 

The  space  and  time  relations  in  equation  (29)  can  be  combined  to  obtain 
other  variables  in  dimensionless  form,  such  as: 


Dp  = BpD*  , 


hp  = Bph*  , etc. 


Substitution  from  equations  (29)  and  (30)  into  equation  (1)  yields: 

♦ ®LD*2£.  + /I\  (T  ) = I (T  ) - 

Yp  at*  ax*  {$)  tT**Jz  5 c 


lar-,2  = v zxJh 


L_  CU*)2  + 3 u*v*' 

|_3x*  D*  ay*  D* 


+ ^EfLfv*  . i %*i!^  + 


3x*  + “ Sxo  - Sjjy  . 
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The  terms  for  gravity,  g,  and  the  rotation  of  the  earth,  f,  are 
not  scaled  because  these  are  not  ordinarily  under  the  control  of  the 
experimenter.  The  terms  for  stress  and  atmospheric  pressure  are  not 
scaled  in  equation  (31)  because  these  have  not  been  clearly  defined  in 
terms  of  the  primary  variables. 

The  terms  in  equation  (31)  for  which  the  variables  have  been  fully 
converted  to  dimensionless  form  contain  the  constant  factors,  a p8p/Yp 
or  glip/ap.  These  factors  can  be  removed  from  all  terms  specified  in 
dimensionless  form  if  equal;  i.e.,  if  the  time  scale,  y,  is  given  by: 


Y 


P 


= 


(32) 


Combining  equations  (29)  and  (32)  show  that  the  transformation  re- 
quired to  convert  the  horizontal  velocity  components  to  dimensionless 
form  is: 


u = '^p  u*  > v = v*  • (33) 

The  transformation  required  for  horizontal  transports  per  unit  width 
as  used  in  equations  (1),  (2),  and  (3)  is: 

Up  = ep  , Vp  - ^V*  . (34) 

Each  term  in  equations  (1),  (2),  and  (3),  as  it  applies  to  the  proto- 
type, can  be  expressed  by  the  product  of  (gBp/a^)  and  the  corresponding 

term  in  dimensionless  form.  The  ratio  of  model  to  prototype  variables 
may  then  be  determined  by  the  ratio  of  the  dimensional  constants.  In  the 
physical  model  for  Masonboro  Inlet  (App.  B) , 


a 300  ’ B = 60  ' 

Thus,  for  the  Masonboro  model,  if 


a = 

’ ap 

, B = 

^ m i/p  1 

Bp  ’ Y a6  38.73 

lV? 

. Vrn 

= gV2 

1 

uP  = 

' vp 

7.746  ’ 

= gae 

1 

Up  * 

' Vp 

464.758  ’ 

= 

a-1  = 

1 

1.549  * 

(35) 


(36) 

(37) 

(38) 
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The  ratio  of  model  boundary  stresses  to  prototype  boundary  stresses 
is  given  by  B2/a.  Assuming  that  the  bottom  stress  can  be  expressed  in 
the  form: 


(xaA  = K^u 


2 + v2' 
Z VZ- 


1/2. 


.(Uz  + Vz) 


21  V2\ 


as  in  equations  (4)  and  (5)  where  u is  scaled  by  equation  (36),  it  is 
found  that 


hn  g 

Kp  = a 


5 for  the  Masonboro  model. 


(39) 


The  stress  coefficient  for  the  bottom  stress  can  be  increased  in  the 
model  to  simulate  the  correct  scaled  bottom  stress  by  increasing  the 
bottom  roughness.  This  technique  cannot  be  applied  at  the  free  surface 
because  the  roughness  of  the  free  surface  is  determined  by  flow  condi- 
tions and  cannot  be  freely  changed  by  the  experimenter.  If  the  free- 
surface  stress  is  modeled,  the  equality  of  stresses  must  be  achieved 
through  the  use  of  a scaling  law  for  airflow  based  on  the  need  for  dupli- 
cating the  surface  stress.  If  the  wind-stress  coefficient  is  a constant, 
the  scaling  law  for  air  velocities  would  be  given  by: 


h = 6 = i 

Vp  /a  3.464 


for  the  Masonboro  model. 


(40) 


where  vm  and  Vp  represent  the  wind  velocity  in  model  and  prototype. 

In  reality,  modeling  the  surface  wind  stress  over  water  is  significantly 
more  difficult  than  equation  (40)  implies,  and  the  surface  stress  is 
generally  neglected  in  hydraulic  models.  An  exceptional  experiment  in 
which  an  effort  was  made  to  model  wind  stress  has  been  reported  by 
Shemdin  (1972). 

The  scaling  laws  for  the  inertial  terms,  (e)  (discussed  in  Sec.  !l), 
in  equations  (1)  and  (2)  are  automatically  accounted  for  by  the  primary 
scaling  relations.  This  constitutes  a major  advantage  of  physical  models. 


The  Coriolis  force  terms,  fV  and  fll,  result  from  the  rotation  of 
the  earth.  The  speed  of  the  earth's  rotation  is  not  under  the  control  of 
the  experimenter  and  so  the  relative  magnitude  of  this  term  in  the  model 
is  less  than  in  the  prototype  by  the  factor  /fi/a  = 1/38.73  for  the  Mason- 
boro model  when  model  and  prototype  are  at  the  same  latitude.  This  is  not 
significant  for  the  flow  inside  Masonboro  Inlet  where  the  Coriolis  force 
is  unimportant,  but  it  can  be  a significant  factor  at  some  locations.  The 
effect  of  a scaled  Coriolis  force  can  be  simulated  in  the  laboratory  by 
building  the  model  on  a turntable  and  rotating  at  an  appropriate  speed. 
Experiments  with  rotating  models  have  been  described  by  Dalvemey  and 
Fontanet  (1959),  Valembois  and  Bonnefille  (1959),  d'Hieres  (1962),  and 
Rumer  and  Robson  (1968).  The  effects  of  the  earth's  rotation  have  been 
neglected  in  the  Masonboro  model  as  in  most  other  physical  hydraulic 
models . 
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In  discussing  the  scaled  pressure  gradients,  it  is  convenient  to 
express  the  pressure  as  the  height  of  a column  of  water  required  to 
express  the  same  pressure.  That  is, 

PJi  = PgH^  , (41) 

where  is  the  height  of  the  water  column. 

The  term,  (g)  (discussed  in  Sec.  II),  described  in  equations  (1)  and 
(2)  may  be  regarded  as  the  product  of  the  fluid  depth  and  the  atmospheric 
pressure  gradient.  Since  the  depth  and  scaling  factor  of  the  fluid  column 
have  been  discussed,  only  the  pressure  gradient  must  be  discussed  here. 

The  prototype  pressure  gradient  is: 


3£  _ dHb 

3 y ™ Sew* 


cipSSx* 


For  an  undistorted  model,  Bp  = a p would  be  independent  of  the  dimensional 
constants,  and  model  and  prototype  pressure  gradients  would  be  identical. 
For  distorted  models. 


/ 1 3Pn\  / 1_  3Pn\~ 1 _ 6_ 
3x  /m\p  3 x / a 


5 for  the  Masonboro  model 


The  atmospheric  pressure  gradient  is  of  minor  importance  in  tidal 
models,  and  the  authors  are  not  aware  of  any  hydraulic  model  studies 
where  an  effort  has  been  made  to  model  the  atmospheric  pressure  gradient. 


Perturbation  and  Viscous  Phenomena  in  Physical  models. 


Perturbation  and  viscous  phenomena  for  physical  models  are  discussed 
in  relation  to  the  mathematical  expressions  which  describe  them.  One  set 
of  terms  due  to  molecular  viscosity  cannot  be  reduced  sufficiently  in  the 
laboratory  to  permit  proper  scaling.  Fortunately,  the  viscous  forces  are 
directly  important  only  in  a very  thin  layer  near  the  boundaries.  The 
direct  effect  of  viscosity  can  always  be  neglected  in  the  prototype.  If 
the  ratio  of  model  depth  to  prototype  depth  is  not  too  small,  the  viscous 
forces  will  also  have  little  direct  influence  on  model  flows.  The  neces- 
sity of  eliminating  the  importance  of  viscosity  in  the  shallow  parts  of 
the  model  basin  is  a prime  reason  for  using  distorted  scale  models.  Even 
in  distorted  models  which  include  regions  of  tidal  flats,  flooded  at  high 
tide  and  exposed  at  low  tide,  there  will  be  a region  of  advancing  water 
at  floodtide  and  a region  of  receding  water  at  ebbtide  in  which  the 
viscous  boundary  layer  may  play  a significant  role.  However,  the  details 
of  the  flow  in  this  region  are  not  of  primary  importance  in  many  inlet 
problems . 


The  turbulent  stresses  are  more  important  than  the  viscous  stresses 
whenever  the  fluid  layer  is  thicker  than  the  viscous  sublayer;  this  will 
include  all  of  the  channels  and  much  of  the  tidal  flat  region.  The 
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structure  of  turbulence  in  the  prototype  is  affected  by  wind  and  surface 
waves  as  well  as  tidal  flow  and  is  not  sufficiently  understood  to  justify 
an  effort  at  modeling  or  precise  simulation.  The  calibration  effort  de- 
scribed in  connection  with  the  discussion  of  surface  stresses  is  assumed 
to  account  for  the  simulation  of  the  principal  effects  of  turbulence  and 
viscosity  on  the  mean  flow. 

The  radiation  stresses  significantly  affect  the  mean  flow  only  in 
regions  where  the  wave  heights  are  changing  significantly  in  space.  The 
greatest  changes  in  wave  height  in  a small  space,  are  those  due  to  refrac- 
tion, shoaling,  reflection,  and  friction  near  the  coast.  These  stresses, 
unlike  the  turbulent  stresses,  are  not  related  by  simple  functions  of  the 
mean  flow.  Hence,  the  effects  of  these  stresses  on  the  mean  flow  can  be 
simulated  independently  from  the  mean  flow.  Sager  and  Seabergh  (1975) 
and  Seabergh  and  Mason  (1975)  showed  that  the  agreement  between  the 
physical  model  of  Masonboro  Inlet  and  the  prototype  can  be  improved  by 
simulating  the  approximate  wave  conditions  for  the  period  when  the  proto- 
type tidal  data  were  collected. 

Although  short-period  surface  waves  cannot  be  modeled  to  scale  in 
distorted  models,  it  is  possible  to  simulate  with  good  approximation,  the 
major  effects  of  surface  waves  on  the  large-scale  flow.  The  surface  waves 
affect  the  large-scale  flow  through  gradients  in  the  radiation  stresses. 

It  is  not  possible  to  obtain  similitude  for  a surface  wave  field  in  a 
distorted  scale  physical  model. 

Keulegan  (personal  communication,  1976)  has  shown  that  linear  wave 
refraction  patterns  can  be  correctly  simulated,  if  the  scaling  factor 
used  for  the  vertical  dimension  in  modeling  the  large-scale  flow  is  used 
for  both  horizontal  and  vertical  dimensions  when  dealing  with  surface 
waves,  and  if  the  effects  of  reflection,  friction,  and  diffraction  can  be 
neglected.  This  study  shows  that  when  these  conditions  are  satisfied, 
even  the  nonlinear  aspects  of  refraction  and  the  gradient  of  the  radia- 
tion stresses  can  be  simulated  in  distorted  scale  models. 

However,  reflection  and  processes  which  depend  explicitly  on  wave- 
length cannot  be  properly  simulated  in  a distorted  model;  e.g.,  the 
effects  of  surface  waves  on  harbor  resonance  cannot  generally  be  in- 
vestigated in  distorted  scale  models. 

To  mathematically  demonstrate  the  modeling  procedures  discussed  above, 
it  is  convenient  to  begin  with  the  equations  for  infinitesimal  waves.  The 
fluid  motions  resulting  from  small-amplitude  gravity  waves,  propagating  in 
the  direction  of  increasing  x on  the  surface  of  a fluid  which  is  other- 
wise at  rest,  are  described  by  the  equations: 


(44) 
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where 

2A  = wave  height,  trough  to  crest 
L = wavelength 

c = phase  velocity  of  waves,  relative  to  the  water 
D = h - Z or  mean  depth  of  the  water 

When  D/L  > 1/2,  tanh  2n  D/L  = 0.99 6,  and  equation  (47)  may  be  stated 
with  sufficient  accuracy  as: 


c 


2 _ &k 

" 2 it  • 


(49) 


When  D/L  is  small,  e.g.,  D/L  < 1/20,  the  function  tanh  2tt  D/L  may 
be  represented  by  the  first  term  in  its  Taylor  series  expansion,  2 it  D/L, 
and  equation  (47)  takes  the  form: 


c2  = gD  . (50) 

Equations  (47)  and  (50)  show  that  the  phase  velocity  of  waves  decrease 
when  propagated  into  water  of  decreasing  depth.  If  a wave  crest  approaches 
the  shore  obliquely,  the  part  of  the  wave  crest  in  shallow  water  will 
travel  slower  than  the  part  in  deep  water  and  the  crest  will  bend  in  a 
direction  more  nearly  parallel  to  the  bottom  contours.  This  wave  refrac- 
tion may  lead  to  an  increase  or  decrease  in  wave  amplitude,  depending  on 
the  bathymetry  and  concomitant  currents. 

Equations  (47)  and  (50)  also  show  that  the  phase  velocity  is  a func- 
tion of  D/L  and  D.  Thus,  the  phase  velocity  in  the  model  will  corres- 
pond, at  every  point,  to  the  phase  velocity  in  the  prototype,  but  only  if 
both  D and  D/L  in  the  model  correspond  at  every  point  to  D and  D/L 
in  the  prototype.  A necessary  and  sufficient  condition  that  D and  D/L 
will  correspond  to  the  prototype  values  in  the  same  manner  at  each  point, 
is  that  L be  scaled  by  the  same  factor  as  D.  In  this  case,  the  ratio 
D/L  is  the  same  in  both  model  and  prototype. 
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The  ratio,  A/L,  given  in  equations  (44)  and  (46)  for  the  velocity 
components  due  to  waves,  will  be  properly  reproduced  in  the  model  only 
if  A and  L are  scaled  in  the  same  way.  It  is  not  surprising  to  find 
that  the  wave  amplitude  must  be  scaled  like  any  other  vertical  dimension, 
but  this  fact  does  provide  a useful  bonus. 

It  is  well  known  that  the  phase  speed  in  deep  water  is  slightly 
greater  than  indicated  by  equation  (49)  by  a factor  which  depends  on 
the  ratio  (A/L)2,  and  that  the  phase  speed  in  shallow  water  is  slightly 
greater  than  indicated  by  equation  (50)  by  a factor  which  depends  on  A/D; 
e.g.,  see  Lamb  (1932;  Secs.  187  and  250).  Dean  (1974)  showed  that,  in 
general,  the  phase  speed  of  waves  is  a function  of  D,  A/D,  and  A/L. 
Thus,  when  D,  A,  and  L are  scaled  by  the  same  factor,  the  phase  speed 
for  waves  of  finite  amplitude  is  correctly  simulated  in  the  distorted  wave 
model.  Chu  (1975)  showed  that  refraction  patterns  can  be  significantly 
altered  by  the  nonlinear  effects  due  to  finite  amplitude  waves. 

The  curvature  of  wave  rays,  for  waves  subject  to  refraction,  is  given 

by: 


39_  _ 1 3c 

3s  c 3n  ’ 


(51) 


where  9 is  the  direction  of  the  wave  ray,  s the  distance  along  the 
wave  ray,  and  n the  distance  along  the  wave  crest.  Since  s and  n 
are  horizontal  directions,  both  are  scaled  in  the  same  manner  and  refrac- 
tion patterns  and  wave  amplitudes  are  properly  reproduced. 

Equation  (50)  shows  that  the  square  of  the  phase  speed  is  scaled  by 
the  vertical  scaling  factor.  Thus,  the  phase  speed  of  the  model  wave  is 
scaled  by  /ff,  which  is  identical  with  the  scaling  factor  for  currents 
(eq.  36).  Therefore,  the  effects  of  currents,  water  depth,  and  wave 
amplitude  on  wave  refraction  can  be  simulated  in  a distorted  scale  model, 
provided  the  vertical  scaling  factor  is  used  for  wavelengths. 

When  D/L  < 1/20,  the  wavelength  has  only  a minor  effect  on  the  phase 
velocity  of  waves  and  it  is  not  necessary  to  use  the  vertical  scaling 
factor  for  scaling  waves  whenever  D/L  in  the  model  is  less  than  1/20. 
For  the  Masonboro  model,  this  corresponds  to  a ratio  of  1/100  in  the 
prototype. 

A predictable  error  may  occur  in  simulating  the  effect  of  refraction 
and  shoaling  on  wave  amplitudes  in  distorted  scale  models  because  of  the 
partial  reflection  of  waves  by  underwater  topographic  features. 

Attempts  to  determine  the  refraction  coefficients  by  measuring  changes 
in  wave  height  along  a wave  ray  may  lead  to  underestimates  of  the  refrac- 
tion coefficients  because  the  increase  in  bottom  steepness  by  the  ratio 
8/a  will  tend  to  increase  the  reflection  of  wave  energy  by  bottom  features 
Steep  slopes  are  known  to  be  more  efficient  reflectors  of  wave  energy  than 
gentle  slopes.  According  to  Kajiura  (1963),  the  reflection  coefficient 
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for  a slope  in  underwater  topography  will  be  the  same  as  that  for  a 
vertical  discontinuity  if  X/L  < 1/20,  when  X is  the  horizontal  length 
of  the  slope,  and  the  reflection  coefficient  is  near  zero  if  X/L  > 1. 

Ippen,  Alam,  and  Bourodimos  (1964)  presented  laboratory  data  which  showed 
that  the  reflection  coefficient  decreases  for  increasing  values  of  the 
ratio  2A/L  for  shallow-water  waves. 

More  wave  energy  is  expected  to  be  reflected  back  to  the  open  sea  in 
the  model  than  in  the  prototype,  but  this  effect  will  be  smaller  in  regions 
of  steep  waves  where  the  wave  amplitude  might  be  expected  to  affect  the 
wave  refraction.  Thus,  refraction  coefficients  determined  by  measuring 
wave  heights  in  the  model  may  be  slightly  less  than  the  refraction  co- 
efficients applicable  to  the  prototype  (the  uncertainty  here  is  due  to 
the  lack  of  knowledge  about  the  effect  of  friction  on  refraction  coeffi- 
cients in  either  model  or  prototype).  In  general,  changes  in  both  wave 
height  and  wave  direction  are  simulated  with  useful  accuracy  in  the 
distorted  scale  tidal  model,  if  both  horizontal  and  vertical  dimensions 
of  the  surface  waves  are  scaled  by  the  vertical  scale  ratio  of  the  tidal 
model . 

The  contribution  of  radiation  stresses  to  the  growth  of  the  mean 
current  for  a wave  traveling  in  the  x-direction  is  given  by  equation  (16) 
in  terms  of  the  motions  due  to  waves.  A first  approximation  to  the  de- 
tailed description  of  wave  motions  is  provided  by  equations  (43)  to  (48) . 
The  component  of  this  contribution  in  any  other  direction  is  given  by  the 
product  of  the  component  in  the  direction  of  wave  travel  and  the  cosine  of 
the  angle  between  the  direction  of  wave  travel  and  the  desired  direction. 

Substitution  from  equations  (29)  and  (43)  to  (48)  into  equation  (16) 
shows  that  the  ratio  of  model  radiation  stresses  to  prototype  radiation 
stresses  is  given  by  82.  Hie  horizontal  gradient  of  the  radiation  stress, 
which  interacts  with  the  mean  flow,  must  involve  the  horizontal  scale  of 
the  mean  flow.  Thus,  the  ratio  of  the  gradient  of  the  radiation  stress 
model  to  prototype  is  given  by  B2/a.  This  is  identical  with  the  ratio 
of  the  surface  stresses,  and  shows  that  with  the  exception  of  effects  due 
to  reflection,  diffraction,  and  friction,  the  effects  of  radiation  stresses 
on  the  mean  flow  are  correctly  simulated  in  distorted  models,  if  the  verti- 
cal scaling  factor  is  used  as  the  horizontal  scaling  factor  for  surface 
waves,  even  though  similitude  is  not  obtained  for  the  surface  wave  field. 

In  many  practical  cases,  the  effects  of  radiation  stress  on  the  mean  flow 
are  secondary  to  the  effects  of  both  friction  and  the  pressure  gradient. 
Most  of  the  effects  of  reflection,  diffraction,  and  friction  are  second- 
ary to  the  effects  of  refraction  in  the  modification  of  surface  waves. 

Thus,  the  effects  of  radiation  stress  on  the  mean  flow  can,  in  most  cases, 
be  simulated  with  good  approximation  in  distorted  scale  models. 

4 . Boundary  Conditions. 

Boundary  conditions  are  easily  established  at  fixed  boundaries  where 
no  flow  occurs.  The  movable  boundary  where  no  flow  crosses  the  boundary 
may  be  established  in  physical  models  in  a natural  way.  TTie  waterline 
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moves  up  and  down  the  beach  and  across  tidal  flats  as  the  water  level 
rises  and  falls.  Complete  similitude  is  not  obtained  because  the  viscous 
boundary  layer  can  occupy  a substantial  part  of  the  thin  water  layers, 
less  than  2.5  centimeters  (1  inch)  in  thickness  above  flat  beaches  and 
tidal  flats  in  the  model  but  not  in  the  prototype.  Surface  tension  will 
cause  the  advancing  waterline  to  develop  a prominent  convex  meniscus  in 
the  model  but  not  in  the  prototype.  These  are  comparatively  minor  short- 
comings of  the  model  and  are  not  expected  to  produce  a significant  dis- 
tortion of  the  flow  away  from  the  affected  areas. 

Riverflow  above  tidewater  can  also  be  reproduced  without  difficulty. 
Tidal  flows  in  rivers  or  channels  leave  the  model  in  an  upstream  direction 
and  create  a greater  problem.  In  the  physical  model  of  Masonboro  Inlet, 
the  effects  of  inland  tidal  flows  are  simulated  by  means  of  reservoirs 
beyond  the  modeled  part  of  the  inlet  that  accept  the  excess  water  during 
floodtide  and  return  it  during  the  ebbtide.  The  volume  and  length  of  the 
reservoir  and  constrictions  on  the  flow  can  be  adjusted  as  part  of  the 
calibration  process  to  improve  the  agreement  of  model  and  prototype  tidal 
heights  and  currents.  Secondary  tide  generators  may  be  used  to  impose 
the  boundary  condition  on  some  inland  tide  channels. 

An  inlet  model  must  also  include  a small  section  of  the  ocean  with 
artificial  boundaries.  The  astronomical  tide  is  simulated  in  physical 
models  by  adding  or  withdrawing  water  through  diffusers  located  at  the 
seabed  near  part  or  all  of  the  artificial  boundaries  in  the  ocean.  The 
waterflow  is  controlled  to  reproduce  the  tidal  level  at  some  designated 
control  point  which  will  likely  produce  the  correct  total  flow  through 
the  inlet.  This  manner  of  controlling  the  flow  through  the  inlet  is 
unlikely  to  have  much  affect  on  the  currents  in  the  inlet  and  in  inland 
channels.  Generally,  the  tidal  currents  in  the  ocean  part  of  the  model 
are  not  properly  modeled  by  this  procedure.  Evidences  of  eddy  formation 
in  tne  ocean  part  of  the  model  should  be  checked  by  prototype  measurements 
before  being  accepted  as  natural  features  of  the  flow. 

IV.  MODELING  TIDAL  INLETS:  NUMERICAL  MODELS 


1.  General . 

Numerical  models  of  tidal  flows  are  constructed  by  representing  the 
elevation  of  the  seabed  and  free  surface  above  a datum  plane,  the  fluid 
transport  vector  or  the  horizontal  velocity  vector,  and  other  variables 
of  primary  interest,  by  their  numerical  values  at  a discrete  set  of  com- 
putation points  at  discrete  times;  i.e.,  the  hydrodynamic  variables  are 
represented  by  numerical  analogs  at  selected  points  in  space  and  time. 

It  is  unnecessary  to  express  the  value  of  each  variable  at  each  point. 

Since  no  physical  flows  are  involved,  there  is  no  need  to  scale  any 
variable.  Accounting  for  the  rotation  of  the  earth,  atmospheric  pressure 
gradients  and  prescribed  wind  stresses  do  not  cause  fundamental  difficulty, 
but  modeling  small-scale  processes  can  be  expensive  because  of  the  great 
spatial  resolution  involved. 
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The  time  resolution  required  for  accurate  calculations  is  nearly 
always  greater  than  needed  for  the  solution  of  engineering  problems,  and 
increases  as  the  resolution  in  space  is  increased.  A reduction  of  the 
space  increment  by  a factor  of  two,  generally  also  means  a reduction  in 
the  time  increment  by  a factor  of  two,  and  for  two  horizontal  dimensions, 
an  increase  in  the  number  of  calculations  required  by  a factor  of  eight. 
Consequently,  small-scale  processes  are  usually  simulated  in  operational 
numerical  models  through  a parameter  relating  the  small-scale  process  to 
a large-scale  variable.  The  representation  of  friction  as  a constant 
multiplied  by  the  square  of  the  velocity  is  an  example  of  this  type  of 
simulation. 

Numerical  models  have  the  distinct  advantage  of  permitting  the  elimi- 
nation of  unwanted  processes,  perhaps  useful  in  determining  the  relative 
importance  of  various  related  phenomena. 

The  ability  to  store  a numerical  model  in  a box  of  punchcards  or  on 
a roll  of  magnetic  tape  when  not  needed  may  also  be  useful.  Since  some 
physical  processes  must  be  simulated  in  numerical  models  for  economy  and 
as  the  simulation  is  not  exact  modeling,  calibration  is  required  for 
numerical  models  as  well  as  physical  models. 

Numerical  models  may  be  constructed  along  many  lines,  and  each  has 
advantages  and  disadvantages.  The  most  common  model  for  two-dimensional 
hydrodynamic  problems  (also  the  one  used  in  both  Masonboro  Inlet  two- 
dimensional  numerical  models)  is  often  called  a finite-difference  model. 
The  region  of  interest  is  covered  with  a rectangular  grid  of  equally 
spaced  points  in  two  orthogonal  directions.  The  computation  points  in 
this  scheme  are  usually  called  "grid  points"  or  "mesh  points."  The  dis- 
tance between  grid  points  is  usually  the  same  value.  As,  in  both  direc- 
tions. Computations  are  generally  made  for  constant  increments  in  time, 
At.  A variety  of  numerical  models  suitable  for  the  study  of  two- 
dimensional  problems  in  hydrodynamics  was  discussed  by  Roache  (1972). 
Reduced  but  constant  values  of  both  As  and  At  may  be  used  in  some 
subregions  of  the  basin  to  provide  greater  detail. 

Natural  basins  are  seldom  represented  accurately  in  terms  of  the 
rectangular  grids  preferred  for  numerical  models.  Curvilinear  boundaries 
are  often  represented  by  line  segments  and  many  right-angle  corners.  Each 
of  these  corners  becomes  the  source  of  false  disturbances  which  may  propa- 
gate throughout  the  system.  In  spite  of  the  disturbances  growing  from 
the  false  corners,  useful  information  has  been  obtained  from  models  using 
"stair-step"  boundaries. 

Jelesnianski  (1976)  and  Wanstrath,  et  al.  (1976)  worked  with  coordina- 
tion transformations  which  permit  some  curvilinear  boundaries  and  unequal 
space  increments  to  be  mapped  onto  rectangles  with  constant  spa!ce  incre- 
ments without  significant  error.  These  transformations  permit  improved 
representation  of  the  boundary  conditions  and  improved  resolution  of 
physical  processes  where  it  is  most  needed. 
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Approximations  of  the  governing  differential  equations  by  the  method 
of  finite  differences  may  be  solved  explicitly  or  implicitly.  Explicit 
models  are  based  on  prediction  equations  which  express  the  future  state 
of  each  variable  explicitly  as  functions  of  present  or  past  values  of  the 
same  and  other  variables.  Implicit  models  are  based  on  the  simultaneous 
solution  of  equations  relating  past,  present,  and  future  states  of  the 
variables.  Each  solution  is  discussed  later  in  this  section. 

An  alternative  model  for  two-dimensional  modeling  of  hydrodynamic 
phenomena,  called  a finite-element  model,  is  gaining  popularity.  For 
many  situations,  the  finite-element  approach  permits  greater  freedom  in 
varying  spatial  resolution  to  satisfy  the  need  for  detail  and  an  accurate 
representation  of  the  boundary  conditions.  The  finite-element  technique 
has  been  discussed  by  Zienkecwicz  (1971) , Connor  (1972) , and  Wang  and 
Connor  (1975).  The  finite-element  technique  has  not  been  used  in  any  of 
the  Masonboro  studies  reviewed  in  this  report.  Hinwood  and  Wallis  (1975a, 
1975b)  provided  a classification  scheme  for  numerical  models  of  tidal 
waters,  along  with  a concise  review  of  141  reports  on  numerical  models 
of  tidal  basins.  Discussions  of  this  report  have  been  provided  by 
Abraham  and  Karelse  (1976)  and  Abbott  (1976). 

2.  Explicit  Prediction  Models. 

When  the  differential  expressions  in  the  motion  equations  (1) , (2) , 
and  (3)  are  represented  by  the  differences  between  the  appropriate  varia- 
bles at  neighboring  grid  points  at  the  same  or  nearby  times,  predictions 
may  be  made  for  a designated  grid  point  "A"  at  time  t + At,  in  terms  of 
conditions  at  the  same  and  neighboring  grid  points  at  times  t and  t - At 
Each  grid  point  may  be  considered  in  turn,  by  systematically  considering 
each  grid  point  in  each  row.  This  procedure,  which  was  used  in  both  two- 
dimensional  models  of  Masonboro  Inlet,  produces  an  explicit  prediction 
model. 

Explicit  prediction  models  have  been  found  to  be  useful  tools  in 
solving  the  two-dimensional  problems  of  meteorology  and  oceanography,  and 
should  be  equally  useful  in  coastal  engineering  if  the  lessons  learned  by 
hydrodynamicists  in  other  disciplines  are  fully  utilized. 

Some  difficulties  and  pitfalls  are  inherent  in  the  numerical  modeling 
art.  One  success  may  lead  to  excessive  confidence  and  misuse  of  the 
technique,  or  one  failure  may  lead  to  distrust  of  a useful  but  misunder- 
stood tool. 

The  best  known  pitfall  in  the  use  of  explicit  prediction  numerical 
models  is  the  development  of  computational  instability  when  At/As  is 
too  large.  Computational  instability,  which  makes  its  presence  known  by 
a growth  in  the  mechanical  energy  as  computed  by  the  model  when  no  corres- 
ponding growth  of  energy  is  expected  in  the  prototype  flow,  may  also 
arise  from  the  nonlinear  terms  in  the  equations  or  from  the  boundary 
conditions . 
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Other  difficulties  are  associated  with  the  development,  in  the  model, 
of  spurious  waves  which  have  no  counterpart  in  the  prototype,  and  a ratio 
of  phase  speed  in  the  numerical  model  to  phase  speed  in  the  prototype 
which  is  less  than  unity.  Both  of  these  difficulties  are  most  serious 
when  the  ratio  of  L/As  is  small,  where  L is  the  wavelength  of  the 
disturbance. 


Round-off  error  which  is  due  to  the  calculations  made  with  a finite 
number  of  significant  digits,  may  become  significant  after  some  large, 
but  finite  number  of  time  steps. 


Serious  growth  of  errors  due  to  any  of  these  causes  can  often  be 
avoided  by  additional  programing,  and  misinterpretation  of  the  errors 
that  occur  can  also  be  avoided  by  understanding  the  nature  and  cause  of 
the  errors.  Therefore,  each  of  these  difficulties  is  discussed  here  in 
detail  to  indicate  probable  trouble  areas  and  the  possible  need  for  cor- 
rective procedures,  not  to  teach  construction  of  perfect  numerical  models. 
The  computational  stability  problem  has  been  discussed  by  Haltner  (1971, 
Chs . 5 and  12)  in  an  introductory  text  in  numerical  weather  prediction 
for  meteorologists. 


The  following  discussion  of  computational  problems  is  designed  to 
permit  the  engineer  to  exercise  some  qualitative  judgment  about  results 
obtained  with  numerical  models  in  general  and  explicit  prediciton  tidal 
hydraulic  models  in  particular. 


a.  Fourier  Series  Representation.  Any  function  which  is  specified 
only  for  2N  discrete  values  of  x can  be  expressed  without  loss  of 
accuracy  as: 
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where  x = nAx,  and  a n and  bn  are  constants  which  can  be  determined 
from  the  tabulated  data  by  standard  techniques.  This  expression  yields 


N + 1 values  of  a^,  and 
b$  are  not  needed  because 
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Equation  (52)  may  be  rewritten  as : 
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Equation  (53)  shows  that  t|>(x)  is  fully  described  by  the  mean  value 
and  harmonics  of  the  interval  2NAx.  Thus  no  perturbation  in  ij^(x)  with 
a characteristic  length  less  than  2Ax  can  be  identified  in  the  available 
data.  Any  perturbation  with  a wavelength  (2  - e)  Ax,  where  c < 2, 
which  does  exist  in  i|>(x) , will  appear  in  the  tabulated  data  to  have  a 
wavelength  of  (2  + e)  Ax  (Fig.  2).  In  the  figure,  wavelengths  shorter 
than  2Ax  appear  in  the  plotted  data  as  having  a longer  wavelength.  This 
is  known  as  a long  wave  alias. 

In  many  two-dimensional  models,  it  is  possible  for  waves  to  be  orga- 
nized along  a diagonal  of  the  computational  grid  squares.  These  waves 
can  be  recognized  in  the  computations  only  by  projections  on  the  compu- 
tational lines.  Consequently,  a wavelength  of  2As  traveling  along  a 
diagonal  will  appear  to  have  a wavelength  of  »^2As,  and  will  be  aliased 
into  a wavelength  of  [2  + (2  - /2) ] As  = 2.586  As.  The  length  used  for 
the  space  averaging,  needed  to  separate  the  tidal  scale  flow  from  pertur- 
bations (discussed  in  Sec.  2 and  App.  B) , must  exceed  2/2As.  However, 
the  critical  length  may  be  different  in  some  models. 

An  early  explanation  of  computational  instability  in  explicit  finite- 
difference  numerical  models,  although  not  entirely  rigorous  but  provides 
much  insight  and  nearly  always  leads  to  the  correct  conclusion,  is  based 
on  the  concept  that  in  any  physical  system  it  is  possible  for  conditions 
at  any  designated  point  A to  be  affected  by  events  at  a second  point 
B,  in  the  time  interval  required  for  energy  to  travel  from  B to  A by 
the  fastest  wave  which  the  physical  system  permits.  This  principle  also 
applies  to  explicit  finite-difference  numerical  computation  models  where 
the  term  "physical  system"  is  replaced  by  "numerical  model."  Correct 
predictions  cannot  be  expected  from  a prediction  equation  which  does  not 
include  all  of  the  information  which  can  affect  the  outcome  of  the  pre- 
diction. This  means  that  the  time  interval.  At,  involved  in  a predic-  "v. 

tion  must  not  exceed  the  time  required  for  the  fastest  wave  to  travel 
from  the  boundary  of  the  influence  region,  to  the  point  for  which  the 
prediction  is  made.  For  physical  models  the  outer  limit  of  the  influence 
region  will  coincide  with  the  circle  of  radius,  r,  centered  on  A, 
where 


r = &tcmax  , 

and  craaa.  is  the  speed  of  the  fastest  wave  which  the  system  can  support. 
For  one-dimensional  numerical  models  the  extreme  limits  of  the  influence 
region  are  given  by: 


x = A ± Atc^  . 

In  numerical  finite-difference  models,  it  is  generally  convenient  to 
use  the  constant  space  increment.  As,  for  both  horizontal  directions. 

This  procedure  was  followed  in  both  of  the  modes  involved  in  this  study. 
Thus,  for  one-dimensional  numerical  models  the  time  increment.  At,  must 
not  exceed  l^s/cJnax.  When  two-dimensional  numerical  models  are  considered, 
it  is  necessary  to  recognize  the  possibility  of  plane  waves  traveling 
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along  a diagonal  of  the  computation  matrix.  In  this  case,  a wave  front 
which  coincides  with  one  diagonal  of  the  computational  square  will  gen- 
erally travel  to  the  nearest  computational  point  on  the  other  diagonal, 
at  a distance  of  (1/2)  /2As  in  a time  interval  At.  For  explicit  two- 
dimensional,  finite-difference  models,  the  maximum  permissible  value  of 
At  is  given  by: 


At  < 
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2c 


max 
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max 


For  tidal  hydraulics  models  involving  incompressible  fluids,  crnax 
generally  given  by: 


cmax  = max  > C55) 

where  D mxx  is  the  maximum  value  of  the  total  depth.  If  the  compressi- 
bility of  water  was  included  in  the  governing  equations,  the  speed  of 
sound  would  have  to  be  used  for  crn2X. 

b.  Implicit  Prediction  Models.  The  implicit  prediction  model  is 
an  alternative  procedure  where  the  prediction  equations  for  all  variables 
at  all  grid  points  are  considered  as  a set  of  N simultaneous  equations 
in  N unknown  values  at  time,  t + At.  With  this  scheme,  all  of  the 
information  available  to  the  model  is  used  in  the  prediction  for  each 
grid  point.  Therefore,  the  calculation  can  be  expected  to  be  uncondi- 
tionally stable  without  regard  to  the  ratio  At/As.  It  can  be  shown 
rigorously  that  this  expectation  is  correct  for  linear  equations  with 
constant  coefficients.  The  implicit  calculation  procedure  results  in  an 
equation  involving  an  N by  N matrix  where  N may  range  from  1,000  to 
10,000.  The  matrix  formed  in  this  way  is  unusually  simple  and  is  gener- 
ally of  the  tridiagonal  or  pentadiagonal  form;  i.e.,  all  elements  more 
than  one  or  two  positions  off  the  main  diagonal  are  zero.  The  solution 
is  usually  obtained  by  iterative  methods.  Sielecki  (1968)  reported  an 
investigation  of  several  iteration  procedures,  none  of  which  led  to  an 
economy  over  the  explicit  prediction  model  used  as  a standard  in  an  in- 
vestigation of  storm  surges. 

Roache  (1972,  pp.  87-95)  discussed  a procedure  which  combined  some 
advantages  of  both  explicit  and  implicit  methods.  In  this  technique, 
called  an  alternating  direction  implicit  (ADI)  model,  the  values  of  h 
are  computed  at  t = nAt,  and  the  values  of  U and  V are  computed  for 
t = (n  + l/2)At,  where  n = 1,  2,  3,  ...  . The  values  of  h and  U for 
future  times,  and  for  all  values  of  x,  are  computed  separately  for  each 
value  of  y.  In  this  calculation,  V is  treated  as  independent  of  time. 
An  implicit  calculation  is  made  for  all  values  of  x between  boundaries 
at  that  value  of  y.  Two  independent  applications  of  the  implicit  inte- 
gration procedure  will  be  made  if,  for  that  value  of  y,  there  is  an 
obstruction  to  the  flow,  as  by  an  island,  with  water  on  both  sides.  When 
all  calculations  for  U are  completed,  the  operation  is  repeated  for  V 
and  h,  assuming  that  U is  independent  of  time,  and  considering  each 
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value  of  x as  one  or  more  independent  problems.  This  procedure  permits 
the  use  of  longer  time  steps  for  a given  space  increment  and  may  be  used 
to  reduce  the  cost  for  a given  space  resolution,  or  to  provide  more  reso- 
lution in  space  at  a fixed  cost.  In  general,  this  procedure  provides  less 
resolution  in  time  and  less  accuracy  for  identical  space  resolution. 

The  loss  in  time  resolution  and  accuracy  is  trivial  in  many  applica- 
tions but  significant  in  others.  No  models  of  this  type  were  included 
in  the  Masonboro  Inlet  comparison  project,  but  Butler  and  Raney  (1976) 
apply  this  procedure  to  the  Masonboro  data.  They  report  acceptable 
results  with  time  steps  30  times  longer  than  would  be  permissible  with 
a purely  explicit  model  and  reduction  in  running  time  by  a factor  of  15. 

Leendertse  (1967,  1970)  has  been  a leading  advocate  of  implicit  and 
hybrid  models  in  the  solution  of  hydrodynamical  problems  of  the  coastal 
region.  He  has  cited  the  unconditional  stability  of  the  solutions, 
demonstrated  for  linearized  equations  as  the  reasons  for  this  preference. 
However,  Weare  (1976)  showed  that  the  nonlinear  terms  used  by  Leendertse 
in  his  numerical  model  but  not  in  his  computational  stability  test,  can 
lead  to  instability  in  many  practical  applications.  The  destabilizing 
effects  of  nonlinear  terms  in  the  prediction  equations  are  discussed 
later.  Leendertse  and  Liu  (1975b)  stated  that  an  explicit  prediction 
scheme  was  used  in  their  latest  model.  Butler  and  Raney  (1976)  acknowl- 
edge the  tendency  for  nonlinear  instability  in  their  model  and  reduce 
the  growth  of  the  instability  by  the  use  of  numerical  filters  which 
selectively  dampen  the  short-period  oscillations  which  make  the  greatest 
contribution  to  the  growth  of  instability. 

Sobey  (1970)  compared  four  finite-difference  schemes  widely  used  in 
the  analysis  of  tidal  flows  in  the  coastal  region  and  showed  that  each 
underestimates  the  phase  speed  of  the  waves.  The  extent  of  the  under- 
estimation decreases  with  decreasing  values  of  the  ratio  As/L,  where 
L is  the  wavelength.  Haltner  (1971,  Ch.  6)  showed  that  the  simple 
numerical  schemes  that  are  most  satisfactory  in  other  respects  generally 
underestimate  the  phase  speed.  Even  the  sign  of  the  phase  speed  may  be 
reversed  in  some  models  when  L = 2As,  but  the  error  is  small  when 
L > 20As.  Schemes  which  give  excessive  phase  speeds  can  be  developed. 
Roache  (1972)  discussed  a scheme  in  which  a correct  mean  phase  speed 
is  maintained  by  alternating  between  procedures  which  underestimate  or 
overestimate  the  phase  speed. 

Haltner  (1971,  Ch.  6)  showed  that  some  simple  computational  schemes, 
which  are  otherwise  satisfactory,  generally  yield  two  wavelike  solutions: 
(a)  The  physical  mode  which  corresponded  to  the  prototype  wave,  and  (b) 
the  computation  mode  which  moved  in  the  opposite  direction  to  the  proto- 
type wave  and  changed  phase  with  each  time  step.  Computational  modes 
are  not  present  in  all  models;  however,  wheg  present,  the  wave  amplitude 
required  to  satisfy  initial  conditions  is  divided  between  the  two  modes. 
The  computational  mode  may  be  dominant  when  L = 2As.  Relative  amplitude 
of  the  computational  mode  decreases  for  longer  waves.  The  computational 
mode  may  be  negligibly  small  for  L > lOAs. 
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If  the  product  of  i|>(x)  and  F (x)  appears  in  the  modeling  equations, 
<Kx)  may  be  expressed  as  in  equation  (53),  and  F(x)  may  be  expressed 
as : 

F(x)  » 2 Bn,  cos  L n_'l  ) . (56) 

n ' = 0 v n ' / 

The  product,  tp  (x)  F(x),  may  be  expressed  as: 

N N , 

*Cx)  F(x)  = 2 2 A^,  cos  , (57) 

n=0  n'  = 0 ' v n'/ 

n,  n'  =0,  1,  2,  3,  . . . N . 

A well-known  trigonometry  identity  may  be  used  to  transform  equation 
(57)  into 


i J>(x)  F (x) 


(58) 


All  terms  of  the  type  (n  - n')  ir/N  may  be  expressed  in  the  form  jir/N, 
where  j < N.  Thus,  each  of  these  terms  corresponds  to  a harmonic  of  the 
computation  region  which  is  adequately  represented  by  the  model.  When 
the  terms  involving  n + n'  are  expressed  in  the  same  manner,  j may 
exceed  N,  if  either  n or  n'  exceeds  N/2.  Equation  (58)  describes 
a flow  of  momentum  from  the  perturbations  of  largest  size  into  perturba- 
tions of  both  greater  and  shorter  wavelengths.  This  is  a natural  process 
which  is  not  well  described  by  the  numerical  models.  Perturbations  with 
wavelengths  greater  than  2As  can  be  properly  represented  on  the  compu- 
tational grid;  perturbations  with  wavelengths  less  than  2As  cannot  be 
properly  represented  and  will  be  aliased,  so  that  they  appear  with  wave- 
lengths longer  than  2As.  The  error  arising  in  this  manner  in  a single 
time  step  is  generally  small,  and  many  time  steps  may  be  required  for  a 
significant  error  to  develop;  however,  errors  from  this  process  may  ulti- 
mately destroy  the  validity  of  the  solution.  This  growth  of  errors  may 
be  reduced  by  the  use  of  digital  filters  designed  to  reduce  or  eliminate 
all  irregularities  with  wavelengths  of  4As  or  less  from  the  data  before 
forming  products.  More  elaborate  filtering  procedures  are  required  for 
long  prediction  times  than  for  short  prediction  times. 


Round-off  error  first  appears  in  the  least  significant  digit  of  the 
numbers  representing  any  of  the  variables.  The  occurrence  of  round-off 
error  is  basically  random.  Thus,  sign  changes  are  frequent,  the  wave- 
lengths and  periods  are  short.  Wavelengths  of  2As  and  periods  of  2At 
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are  most  common.  Round-off  errors  grow  if  the  computation  scheme  is  in 
any  way  unstable.  The  growth  may  be  exponential  with  the  number  of  time 
steps.  With  modern  computers  having  a word  length  of  32  bits  or  more, 
several  hundred,  or  sometimes  several  thousand  time  steps  are  usually 
required  before  round-off  error  becomes  significant.  Round-off  error 
can  generally  be  reduced  by  increasing  the  word  length  used  in  the  cal- 
culations. The  growth  of  round-off  error  will  ultimately  overcome  the 
effect  of  a straightforward  reduction  of  the  space  increment.  As,  in 
improving  the  accuracy  of  the  calculations,  especially  when  minicomputers 
with  short  word  lengths  are  used.  However,  large  numbers  of  calculations 
are  generally  cheaper  on  modem  computers  with  long  computational  word 
lengths,  so  the  difficulty  here  may  be  more  a matter  of  convenience  than 
economics  and  can  often  be  solved  by  using  a computer  with  a longer  compu- 
tational word.  The  total  cost  of  the  computation  will  increase  when  in- 
creased accuracy  is  obtained  solely  by  a reduction  in  the  grid  increment. 


Finite-difference  equations  are  formed  by  using  finite  differences 
to  approximate  the  differentials  which  appear  in  the  governing  differen- 
tial equations.  The  error  resulting  from  this  approximation  is  called 
the  truncation  error,  and  can  be  estimated  as  shown  below.  In  forming 
an  expression  for  3y/3x,  at  x = nAx,  y is  expanded  as  a Taylor 
series  around  nAx.  Thus, 

yfnAx  + Ax)  = y(nAx)  + y' (nAx)  Ax  + y"(nAx)  Ax2  + ...  . 


By  rearranging  terms , 


h. 

3x 


nAx 


y(nAx  + Ax)  - y(nAx) 

Ax  1 ’ 


where 


1 ifx 

2 3x2 


Ax  + 
nAx 


The  above  equation  is  of  first  order  because  it  depends  on  the  first 
power  of  Ax.  An  approximation  which  is  accurate  to  second  order,  e.g., 
the  error  depends  on  the  second  power  of  Ax,  may  be  formed  by  combining 
expansions  for  y(nx  ± Ax)  to  obtain: 


where 


3J1 

3x 


nAx 


y(nAx  + Ax)  - y(nAx  - Ax) 
2Ax 


e2  » 


L*h. 

3 3x3 


3x2  . 
nAx 
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Higher  ordered  approximations  may  be  derived  by  considering  more  values 
of  y and  eliminating  more  derivatives  from  the  resulting  expansions. 
By  assuming  that  y has  the  form. 


y = cos 


higher  ordered  approximations  can  be  shown  to  give  better  estimates  of 
the  derivatives  when  the  wavelength  is  long  relative  to  the  interval 
used  in  defining  the  derivatives.  However,  for  the  shortest  waves  sup- 
ported by  the  system,  the  use  of  higher  ordered  approximations  may  lead 
to  less  accurate  estimations  of  the  derivatives  than  low-ordered  approxi- 
mations . 


c.  Error  Control.  Numerical  models  contain  many  sources  of  errors. 
Nearly  all  originate  with  small-scale  perturbations  that  cannot  be  re- 
solved on  the  computational  grid,  or  with  the  smallest  perturbations  that 
can  occur  in  the  grid.  The  probability  of  significant  errors  from  any 
source,  other  than  round  off,  can  be  reduced  by  reducing  the  space  incre- 
ment, As,  to  obtain  a better  resolution  of  the  phenomena  being  studied. 
Increasing  the  resolution  will  increase  the  number  of  time  steps  required, 
thus  increasing  the  cost  of  the  calculation  and  increasing  the  probability 
that  round-off  error  will  become  significant. 

The  probability  of  significant  error  may  also  be  reduced  by  removing 
the  small-scale  phenomena  from  the  governing  equations  and  the  boundary 
and  initial  conditions  at  the  beginning  of  the  calculation.  Since  the 
product  of  two  variables  may  contain  perturbations  of  smaller  scale  than 
either  variable,  the  stability  of  the  calculations  can  often  be  improved 
through  the  use  of  smoothing  or  filtering  operations  at  various  times 
during  the  calculations.  Even  when  errors  cannot  be  entirely  eliminated, 
it  may  be  helpful  to  know  that  the  large-scale  features  are  more  nearly 
correct  than  small-scale  perturbations.  Small-scale  errors  due  to  faulty 
boundary  conditions  are  often  no  more  than  2 or  3 As  in  length;  there- 
fore, since  the  phase  speed  in  the  numerical  model  is  an  increasing  func- 
tion of  As/L,  the  errors  only  slowly  propagate  into  the  interior  of  the 
computation  region.  Thus,  the  solution  tends  to  be  more  reliable  near 
the  center  of  the  computation  region,  and  least  reliable  near  boundaries 
where  the  true  boundary  condition  may  be  unknown,  or  not  properly  repre- 
sented on  the  computational  grid.  The  Masonboro  Inlet  numerical  models 
have  been  used  only  for  a study  of  tide  propagation  into  the  inlet-bay 
system.  The  length  of  the  tide  wave  is  many  times  the  length  of  the  basin 
considered. 

Initially,  the  flow  variables  contained  no  short  wave  perturbations. 
The  phase  speed  of  tide  wave  in  the  model  should  have  been  within  1 to  2 
percent  of  the  prototype  speed.  Although  the  computational  mode  does  not 
appear  in  any  of  the  models  used  in  this  study,  the  success  of  the  model 
in  reproducing  tidal  flow  should  not  lead  to  an  assumption,  without  fur- 
ther examination,  that  the  same  models  will  provide  satisfactory  solutions 
for  a study  of  harbor  oscillations,  where  prominent  oscillations  with 
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short  wavelengths  when  compared  to  the  entire  basin  may  occur,  although 
both  phenomena  are  governed  by  the  same  differential  equations.  With 
extended  operation  of  the  tidal  cycle  in  the  Masonboro  model,  the  solu- 
tion can  be  expected  to  deteriorate  due  to  nonlinear  interaction  of  the 
tide  with  small  features  of  the  basin.  The  model  should  perform  better 
in  bays  with  more  regular  coastlines  and  which  are  larger  relative  to 
the  space  step  employed.  Thus,  the  model  may  perform  better  when  applied 
to  Galveston  Bay,  Apalachicola  Bay  or  the  North  Carolina  sounds  than  when 
applied  to  Masonboro  Inlet. 

d.  Selection  of  Terms  for  Numerical  Modeling  and  Representation. 

Even  if  the  numerical  modeling  was  accomplished  in  an  ideal  manner,  the 
model  could  give  incorrect  results  if  important  physical  processes  are 
omitted,  or  improperly  represented  in  finite-difference  form.  The  im- 
portance of  various  physical  processes  was  discussed  in  Section  II. 
Difficulties  in  modeling  or  simulating  these  processes  in  numerical 
models  are  discussed  below. 

Nonlinear  numerical  models  always  require  more  computer  time  than 
linear  models  with  the  same  resolution,  although  in  some  cases  the  in- 
crease may  be  small.  Adding  nonlinear  terms  to  a linear  numerical  model 
always  increases  the  probability  of  computational  instability.  There  is 
a distinct  advantage  in  linearizing  the  governing  equations  unless  a 
significant  reduction  in  error  can  be  expected  to  result  from  including 
the  nonlinear  terms.  Specification  of  the  conditions  and  procedures 
which  justify  linearization  in  numerical  models  is  an  essential  purpose 
of  this  discussion. 

The  first  term  in  equation  (1),  3U/9t,  is  clearly  linear  and  re- 

quires no  additional  discussion. 

The  second  term  in  equation  (1)  may  be  rewritten  as: 

g°£-  g(h,  - Z)[l  ♦ (59) 


where  hi  is  the  mean  sea  level  (MSL)  or  mean  water  level  over  an  inter- 
val much  longer  than  any  prominent  period  of  the  phenomena  being  modeled. 
The  product, 
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is  quadratic  in  h.  If  (h  - hi)  « (hi  - Z)  everywhere,  neglecting  the 
quadratic  term  will  not  lead  to  serious  error.  Neglecting  this  term  may 
be  desirable  if  (h  - hi)  < (hi  - Z)  everywhere  and  much  less  at  all  but 
a few  isolated  points,  or  much  less  everywhere  but  the  line  of  grid  points 
nearest  the  shore.  The  quadratic  term  in  h is  sometimes  neglected  in 
the  study  of  tidal  flow  within  fixed  boundaries.  Neglecting  the  quadratic 
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term  when  (h  - hi)  and  (hi  - Z)  are  of  the  same  order  of  magnitude 
over  any  substantial  area  of  the  basin  may  lead  to  significant  error. 

The  quadratic  term  cannot  be  neglected  in  models  where  boundaries 
move  inland  with  rising  water.  This  nonlinear  term  is  included  in  both 
Masonboro  models  because  it  is  necessary  to  account  for  the  flooding  of 
tidal  flats  in  the  model. 

The  bottom  stress,  term  (c)  in  Section  II,  must  be  explicitly  stated 
as  input  information  from  an  external  source  or  simulated  as  a clearly 
defined  function  of  the  physical  variables  in  numerical  models.  The 
technique  of  implicit  simulation  by  adding  artificial  roughness  to  the 
boundaries  used  with  physical  models  is  not  available. 

The  usual  procedure  is  to  express  the  bottom  stress  as: 

(Ts^Z  = K|q|2  cos  0 , (61) 

where  q is  the  magnitude  of  the  vertically  averaged  velocity,  0 the 
angle  between  the  velocity  and  the  x-axis,  and  K a coefficient  assumed 
to  depend  on  bottom  roughness.  The  bottom  stress  is  often  expressed  in 
the  algebraically  equivalent  form  of  equation  (4).  In  early  models,  K 
was  taken  as  a constant  throughout  the  model.  This  saved  computer  space 
and  was  a convenient  choice  while  developing  the  more  important  aspects 
of  the  model.  In  general,  the  calculations  are  not  very  sensitive  to 
the  value  of  K.  In  many  other  models  developed  for  modem  computers 
where  storage  space  is  not  a premium,  K has  been  treated  as  a function 
of  position;  i.e.,  a different  value  of  K may  be  used  for  each  grid 
point.  The  coefficient,  K,  may  be  taken  initially  as  a decreasing 
function  of  depth  and  later  modified  on  a point  by  point  basis  as  needed 
for  model  calibration.  Masch,  Brandes,  and  Reagan  (1977)  follow  this 
procedure  in  a model  of  Masonboro  Inlet. 

Reid  (1957)  analyzed  channel  flows  subject  to  wind  action,  and  showed 
that  the  value  of  K when  used  as  in  equation  (4),  is  a function  of  wind 
stress  as  well  as  bottom  roughness.  Harris  (1967)  reported  that  Reid's 
formulation  of  the  bottom  friction  term  gives  better  predictions  for 
hurricane  storm  surges  on  the  open  coast  than  could  be  obtained  with  the 
form  of  equation  (61).  This  formulation  is  nonlinear  and  tends  to  pro- 
duce a growth  of  small-scale  perturbations  with  the  associated  aliasing 
problems  if  the  running  time  for  the  model  is  large,  unless  eliminated 
by  the  application  of  a smoothing  operation.  Many  investigators,  e.g., 
Dronkers  (1964),  introduced  heuristically  derived  linear  friction  laws 
when  working  with  equations  that  were  otherwise  linear.  However,  the 
quadratic  friction  law  does  not  have  a rigorous  justification,  (see 
Sec.  II). 

Jelesnianski  (1967)  expanded  on  a development  by  Platzman  (1963)  to 
show  that  a linear  friction  law  could  be  derived  for  geophysical  flows  in 
a more  rigorous  manner  than  that  used  in  deriving  the  more  widely  used 
quadratic  law,  and  that  the  application  of  this  law  permitted  a signifi- 
cant improvement  in  hurricane  storm  surge  predictions  along  tne  open 
coast  over  those  obtained  with  Reid's  (1957)  friction  law. 
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Jelesnianski,  in  developing  his  friction  law,  recognized  that  bottom 
friction  is  a result  of  a vertical  diffusion  of  momentum,  and  that  dif- 
fusion requires  finite  time,  more  time  being  required  for  greater  depths. 
Because  of  the  earth's  rotation  during  the  diffusion  time,  the  stress 
vector  at  the  bottom  is  not  generally  parallel  to  the  mean  flow.  How- 
ever, for  short-time  periods  and  shallow  water,  the  bottom-stress  vector 
is  nearly  parallel  to  the  mean  flow.  In  general,  the  velocity  does  not 
vanish  outside  the  viscous  boundary  layer  near  the  bottom. 


The  bottom  stress  in  Jelesnianski ' s model  is  expressed  as: 
v3(u  + iv) 


3z 


z = Z 


s£(u  + iv) 
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z = Z 


where 

v = effective  viscosity  coefficient 

u,  v = velocity  components 

i = >^T 

S£  = slip  coefficient 

Thus,  the  Jelesnianski  friction  law  is  linear  in  the  flow  values  but  re- 
quires two  coefficients.  Jelesnianski  (1967)  reported  that  the  solution 
is  less  sensitive  to  the  exact  values  of  the  coefficients  v and  s, 
than  it  is  to  the  exact  value  of  K when  equation  (4)  or  (61)  is  used 
to  simulate  the  bottom  stress  term.  Forristall  (1974)  extended  the 
analytical  path  of  Platzman  (1963)  and  Jelesnianski  (1967)  to  provide  a 
model  for  computing  the  vertical  profile  of  the  horizontal  currents  at 
selected  locations  without  the  expense  of  using  a three-dimensional  grid. 

The  stress  at  the  free  surface  is  normally  attributed  to  the  wind  and 
must  be  supplied  explicitly  to  the  computer  routine  used  for  the  calcula- 
tions of  water  motions.  The  application  of  an  arbitrary  wind-stress  field 
to  the  numerical  model  does  not  cause  any  difficulty  unless  the  wind-stress 
field  contains  small-scale  features  that  cannot  be  properly  represented  at 
the  grid  scale  used.  Wind  stress  was  considered  by  Masch,  Brandes,  and 
Reagan  (1977)  in  the  model  of  Masonboro  Inlet.  The  ability  to  impose  a 
wind-stress  field  on  the  equations  of  motion  for  waterflow  is  a prime 
reason  for  adopting  numerical  models  for  the  study  of  coastal  phenomena. 

Wind  stress  is  a secondary  aspect  of  most  inlet  problems  and  a con- 
sideration of  wind-stress  effects  was  not  part  of  the  specifications  for 
the  Masonboro  Inlet  study.  Therefore,  the  problems  of  developing  the 
wind-stress  field  from  the  available  meteorological  data  are  not  discussed 
in  this  study. 
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In  physical  models  the  advective  terms,  (e)  of  Section  II, 

3_  U 3_  UV 
3x  D + 3y  D ' 

are  always  treated  implicitly  without  any  special  attention  by  the  experi- 
menter. However,  the  terms  appear  in  numerical  models  only  if  they  are 
explicitly  expressed.  These  nonlinear  terms  may  be  improperly  represented 
in  numerical  models  unless  an  effort  is  made  to  eliminate  the  generation 
of  perturbations  too  small  for  proper  representation  on  the  calculation 
grid.  The  physical  importance  of  these  terms  was  discussed  in  Section  II, 
and  must  be  included  if  omission  is  likely  to  produce  a significant  error. 
The  terms  were  included  in  the  most  successful  of  the  Masonboro  models 
(Masch,  Brandes,  and  Reagan,  1977),  but  were  expressed  approximately  as: 

U 3U  V 3U 
D 3x  + D 3y  * 

Butler  and  Raney  (1976)  and  others  also  used  these  incomplete  forms,  but 
the  advective  terms  were  omitted  by  Chen  and  Hembree  (1977).  The  repre- 
sentation of  the  advective  terms  in  the  form  used  by  Masch,  Brandes,  and 
Reagan  would  represent  an  improvement  over  their  omission  in  some  but  not 
all  practical  conditions.  This  problem  requires  more  study  for  a defini- 
tive evaluation. 

Three  effects  of  the  earth's  rotation  (Coriolis  acceleration)  on 
tidal  flows  were  discussed  under  term  (e)  in  Section  II,  3.  The  most 
important  effect  in  modeling  flows  of  tidal  period  is  the  continual  ac- 
celeration of  the  current  to  the  right  of  the  existing  current  direction 
in  the  Northern  Hemisphere  and  to  the  left  in  the  Southern  Hemisphere 
which  can  be  reproduced  in  numerical  models  without  great  difficulty. 

This  feature  is  included  in  all  large-scale  numerical  weather  models. 

The  effect  of  Coriolis  acceleration  in  producing  a rotation  of  the  mean 
current  with  increasing  depth  and  the  effect  of  Coriolis  acceleration 
turbulent  stress  have  been  successfully  modeled  by  Platzman  (1963), 
Jelesnianski  (1967),  and  Forristall  (1974).  The  authors  are  not  aware 
of  any  experiments  in  which  the  generation  of  the  inertial  circle  was 
simulated  as  an  incidental  feature  of  the  calculation  of  tidal  or  storm 
surge  flows.  It  is  assumed  that  effects  of  the  inertial  circle  on  the 
mean  flow  cannot  be  simulated  by  means  of  parameters  and  mean  flow  vari- 
ables unless  the  averaging  period  extends  over  many  days. 

The  possibility  of  including  the  atmospheric  pressure  gradient  effect 
(term  (g)  in  Sec.  II,  3)  in  a straightforward  manner  is  a minor,  but  some- 
times significant  advantage  of  numerical  models  of  water  motions  in  the 
coastal  zone.  When  atmospheric  pressure  is  considered,  the  essential 
information  must  be  explicitly  supplied  as  input  data.  Hie  atmospheric 
pressure  gradient  acts  as  a correction  to  the  surface  slope  term.  When 
pressure  information  is  available,  it  may  be  treated  in  the  same  manner 
as  the  slope  term  with  regard  to  linearization. 
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Each  term  in  equation  (2)  is  analogous  to  a corresponding  term  in 
equation  (1)  and  all  terms  in  equation  (3)  are  well  specified  and  linear. 
Thus,  no  additional  discussion  of  the  principal  terms  is  required.  In 
numerical  models,  as  in  physical  models,  it  is  impossible  to  treat  the 
perturbation  processes  explicitly;  however,  simulation  of  the  most  impor- 
tant effects  of  the  small-scale  processes  on  the  large-scale  flow  is  often 
practical . 

Turbulent  stresses  act  to  smooth  irregularities  and  eliminate  sharp 
gradients  in  the  prototype  flow.  Smoothing  operators  accomplish  the  same 
effect  in  numerical  models.  In  many  cases,  the  large-scale  flow  is  not 
highly  sensitive  to  the  amount  of  smoothing  involved.  Occasionally,  it 
is  necessary  to  increase  the  amount  of  smoothing  in  some  part  of  the  model 
to  obtain  a good  reproduction  of  the  flow.  Overland  and  Myers  (1976)  re- 
ported on  a one-dimensional  model  of  the  Cape  Fear  River,  North  Carolina, 
in  which  the  horizontal  current  over  the  tidal  flats  was  related  to  the 
flow  in  the  main  channel  through  an  eddy  viscosity  term  in  a form  which 
is  equivalent  to  writing  equation  (11)  as: 

„ . 8u 

s xt  = Az/arjty  • (63) 

Leendertse  and  Liu  (1975a)  described  a three-dimensional  numerical 
model  in  which  the  stress  terms  are  represented  by  horizontal  and  vertical 
eddy-stress  coefficients  of  the  type  described  by  equation  (13). 

The  radiation  stresses  are  generally  considered  in  investigations  of 
wave-generated  currents.  Numerical  models  which  combine  calculations  of 
wave  refraction  and  radiation  stresses  for  the  calculation  of  nearshore 
currents  have  been  reported  by  several  investigators,  e.g.,  Noda,  et  al. 
(1974).  Liu  and  Mei  (1975)  considered  refraction  diffraction  and  radia- 
tion stresses  near  a breakwater.  No  numerical  model  is  known  in  which 
the  effect  of  waves  on  tidal  currents  in  an  estuary  has  been  considered. 

e.  Boundary  Conditions.  Boundary  conditions  are  easily  applied  for 
fixed  boundaries  with  zero  flow  normal  to  the  boundary  in  both  physical 
and  numerical  models.  Movable  boundaries  are  more  difficult  to  treat  in 
numerical  models  than  in  physical  models.  Most  aspects  of  the  numerical 
modeling  of  geophysical  hydrodynamic  phenomena  have  been  extensively 
treated  by  meteorologists  interested  in  numerical  weather  prediction, 
or  other  scientists  interested  in  basic  research.  Numerical  modeling 
of  boundaries  which  change  in  location  due  to  tidal  fluctuations  is  an 
exception.  The  only  treatment  used  in  the  Masonboro  model  studies  as 
part  of  the  General  Investigation  of  Tidal  Inlets  (GITI)  is  the  procedure 
reported  by  Reid  and  Bodine  (1968).  Modifications  of  this  procedure  were 
used  in  both  two-dimensional  models  of  Masonboro  Inlet.  The  boundary 
location  in  the  numerical  model  is  allowed  to  change  discontinuously  by 
the  distance  As  in  either  direction,  as  the  water  rises  above  or  falls 
below  designated  critical  values  which  depend  on  the  elevation  of  the 
land  surface  surrounding  a grid  point.  No  theoretical  evaluations  of 
this  procedure  are  known.  The  scheme  is  described  by  Masch,  Brandes, 
and  Reagan  (1977). 
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Riverflows  above  tidewater  cause  no  great  difficulties  in  either 
physical  or  numerical  models.  Tidal  flows  in  rivers  or  channels  leaving 
the  model  in  an  upstream  direction  cause  difficulties  in  both  numerical 
and  physical  models.  The  method  of  treating  these  inland  tidal  boundaries 
adopted  for  the  physical  model  could  have  been  applied  in  the  numerical 
models,  but  this  would  have  required  detailed  calculation  of  the  water 
motions  in  the  fictional  reservoirs  and  in  the  tidal  channels  and  would 
have  been  more  expensive  for  the  numerical  than  for  the  physical  models. 
Masch,  Brandes,  and  Reagan  (1977)  simulate  tidal  records  at  the  internal 
water  boundaries  in  the  form  of  amplitude  and  phase  changes  applied  to 
the  open  ocean  tide.  Chen  and  Hembree  (1977)  simulate  these  internal 
boundaries  by  extrapolation  from  the  interior  of  the  computation  region. 
The  relative  merits  of  the  two  basic  expressions  which  can  be  used  for 
open  water  boundaries  were  discussed  by  Reid  (1975)  and  Forristall  (1975). 
It  appears  that  no  available  procedure  for  dealing  with  open  water  bound- 
aries in  numerical  models  is  ideal;  however,  the  procedure  used  by  Masch, 
Brandes,  and  Reagan  (1977)  may  be  preferable  to  that  used  by  Chen  and 
Hembree  (1977). 

Both  numerical  and  physical  models  of  inlets  require  artificial  bound- 
aries around  a small  part  of  the  ocean.  In  both  models  the  information 
required  for  controlling  the  flow  is  provided  by  specifying  the  tidal 
height,  as  a function  of  time,  at  a specified  point.  Currents  in  the 
numerical  models  are  controlled  by  the  slope  of  the  water  surface  near 
the  artificial  boundary.  Masch,  Brandes,  and  Reagan  (1977)  and  Chen  and 
Hembree  (1977)  use  different  expressions  for  this  purpose.  The  expres- 
sions by  Reid  (1975)  and  Forristall  (1975)  applied  to  the  tidal  channels 
and  to  the  open  water  conditions  in  the  ocean. 

Controlling  the  currents  by  the  computed  slope  of  the  water  surface 
in  numerical  models  is  less  restrictive  than  the  procedure  used  for 
physical  models,  since  tidal  currents  parallel  to  the  shore  can  develop. 
Cotidal  charts,  which  give  the  relative  phase  of  the  tide  over  an  entire 
ocean  basin,  generally  show  that  the  tide  wave  travels  parallel  to  the 
shore,  indicating  that  tidal  currents  should  have  a component  parallel 
to  the  shore.  However,  the  cotidal  chart  by  Defant  (1961)  showed  that 
much  of  the  U.S.  Atlantic  coast  is  an  exception,  and  that  tidal  currents 
parallel  to  the  shore  are  not  expected  near  Masonboro.  Thus,  the  possi- 
bility of  predicting  tidal  currents  parallel  to  the  shore  is  not  as  im- 
portant for  Masonboro  Inlet  as  for  other  inlets,  or  as  important  for 
tidal  problems  as  for  storm  surge  models. 

V.  A LUMPED  PARAMETER  MODEL 


1.  General . 

The  design  engineer  or  researcher  frequently  needs  a simple  model  to 
get  a "first  look"  at  inlet  hydraulics.  The  model  should  be  easy  to  use, 
inexpensive,  require  a minimum  number  of  input  parameters,  and  give 
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reasonable  results  quickly.  As  in  other  mathematical  or  physical  models, 
this  model  is  optimized  by  adjusting  the  input  data  to  simulate  known 
conditions . 

A mathematical  model  for  tidal  inlets,  designed  to  meet  these  criteria, 
is  the  "lumped  parameter  inlet  model"  proposed  by  Keulegan  (1967).  In  the 
lumped  parameter  model,  all  coefficients  which  could  be  adjusted  to  im- 
prove model  performance  were  combined  into  a single  parameter. 

The  original  development  was  based  on  the  consideration  of  a minimum 
number  of  terms  which  might  possibly  provide  useful  results  for  the  sim- 
plest cases,  and  a recognization  that  the  bottom  roughness  could  not  be 
determined  accurately  without  a detailed  investigation.  Manning's  n 
may  be  uncertain  by  a factor  of  two.  Thus,  many  secondary  terms  may  be 
neglected  without  a significant  increase  in  the  uncertainty  of  the  cal- 
culations . 

The  Keulegan  model  has  led  to  satisfactory  resolution  of  many  inlet 
problems.  In  other  applications,  the  results  ohtained  with  the  Keulegan 
model  have  been  unsatisfactory.  The  model  has  been  extended  by  Huval  and 
Wintergerst  (1977)  and  others  by  adding  new  terms  to  permit  extension  of 
the  model  to  general  conditions. 

The  equations  by  Keulegan  to  describe  the  model  and  the  definition 
sketch  (Fig.  1 in  Keulegan,  1967)  are  inconsistent. 

Physical  insight  for  the  dominant  control  of  tidal  flows  in  many  in- 
lets may  be  gained  by  a revision  of  the  definition  sketch  shown  in  Figure 
3(a). 

The  original  figure  in  Keulegan  (1967)  identified  only  two  points  in 
a horizontal  plane:  the  boundary  between  the  sea  and  the  inlet,  x2,  and 
the  boundary  between  the  bay  and  the  inlet,  x1.  The  figure  and  equations 
by  Keulegan  seem  to  imply  that  the  pressure  gradient  from  x2  to  Xj  is 
responsible  for  the  acceleration  required  to  produce  the  inflow  velocity 
at  x2  during  floodtide,  and  that  the  pressure  gradient  from  x^  to  x2 
is  responsible  for  the  acceleration  required  to  produce  the  inflow  veloc- 
ity at  X!  during  ebbtide.  However,  this  is  not  physically  possible. 

If  two  additional  positions,  xq  and  x3,  are  postulated  (Fig.  3,b) , 
a lower  slope  from  x3  to  x2  (Keulegan 's  AHj)  can  produce  the  accelera- 
tion required  to  generate  the  flood  currents,  and  the  slope  between  x2 
and  Xj  (Keulegan' s AH2)  could  exactly  balance  the  friction  and  kinetic 
energy  terms  as  he  has  postulated.  To  apply  the  same  logic  on  ebbtide, 
it  is  necessary  to  assume  that  the  acceleration  required  for  the  ebb 
current  is  gererated  by  the  slope  from  x0  to  Xj  (Fig.  3,c).  This 
approach  is  excellent  for  presenting  the  concept  that  in  many  inlets, 
the  slope  through  the  inlet  is  approximately  balanced  by  friction  within 
the  inlet,  and  in  many  situations  will  permit  a simple  and  useful  first 
approximation  to  inlet  flows. 
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Figure  3.  Gradient  of  water  surface  over  channel:  (a)  as  presented  by 
Keulegan  (1967),  (b)  as  modified  floodflow  conditions,  and 
(c)  as  modified  ebb  flow  conditions. 
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However,  this  interpretation  does  not  provide  a satisfactory  approach 
to  generating  improvements  to  the  first  approximation,  since  it  does  not 
permit  a straightforward  application  of  the  equation  for  the  conservation 
of  energy,  momentum,  and  mass. 


If  the  purpose  of  the  model  is  quantitative  calculations,  the  figure 
appears  to  be  correctly  drawn  and  the  equations  appear  to  be  missing 
several  terms  which  are  negligibly  small  in  many  practical  applications, 
but  important  in  others.  Most  of  the  material  in  this  section  was 
developed  on  this  assumption. 


A more  rigorous  development  than  that  originally  employed,  which 
clearly  shows  the  implicit  assumptions  used  by  Keulegan  (1967)  and  by 
Huval  and  Wintergerst  (1977)  may  be  developed  by  integrating  the  one- 
dimensional equations  of  motion  from  the  open  sea  to  the  bay  or  lagoon. 
This  approach  discloses  the  terms  omitted  in  the  original  equations.  In 
many  applications  these  terms  do  not  have  a significant  effect  on  the 
calculated  results,  but  are  important  in  others.  The  one-dimensional 
equations  of  motion  are  developed  in  Appendix  C by  integraing  the  two- 
dimensional  equations  of  motion  with  respect  to  y from  one  side  of  a 
channel,  yi(x),  to  the  other  side  of  the  channel,  y2(x). 


The  resulting  equations  (B13)  and  (B14)  in  Appendix  B may  be  stated 


as : 


3u  1 3 =2  3h  1 [ , . , 1 

3t  + 2 3^  + g3l  + Ac  J dy  = AC 


= 2 

u I (R  + E)  dy 


y 2 


|— D + 
3t 


3 K=  . -1  ( r n „ . , D^b  =3b"l 
^Ou  = b J (R  + E)  dy  - + u^J 

yi 


, (64) 
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where 


YZ 


,(x,t)  = J D dy 
yi 


is  the  cross-sectional  area  of  the  inlet  at  position  x along  the  channel 
at  time,  t, 

YZ 

u(x,t)  = A -al  j U dy 

yi 
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Yz  h 

= A^1  j j u dz  dy  . 

yi  z 


Thus,  u(x,t)  is  the  average  of  u(x,y,z,t)  in  that  cross  section, 

Yl 

Rc  = " lx  / [D(u") 2 + 2uD"u"  + D"(u")2]  dy 

yi 


u"(x,y, z,t) 


o-a." 

dX 


yi 

u(x,y,z,t)  - 


u(x,t)  , 


and 


y 2 

h(x)  = b"1  J h(x,y)  dy  , 

y\ 


h' (x,t)  = h - h , 

D'(x,y,t)  = D - 5 , 

b(x.)  = y2  - y1  width  of  channel  . 


Rc  contains  all  terms  which  appear  when  the  one-dimensional  momentum 
equation  is  derived,  without  approximation,  from  the  three-dimensional 
equations,  but  which  are  generally  neglected  in  the  one-dimensional 
momentum  equation. 

The  right-hand  sides  of  equations  (64)  and  (65)  are  closure  terms, 
normally  neglected  in  dealing  with  the  one-dimensional  equations  of 
motion.  These  terms  are  generally  small  and  are  given  primarily  for 
completeness,  but  several  of  the  terms  can  be  locally  significant.  The 
physical  processes  which  the  terms  represent  were  discussed  in  Section  II. 
None  of  the  terms  were  considered  in  the  lumped  parameter  models  of 
Keulegan  (1967)  and  Huval  and  Wintergerst  (1977). 

In  most  applications,  the  displacement  of  the  free  surface  is  assumed 
infinitesimal  with  respect  to  the  total  depth,  and  the  elevation  of  the 
channel  bottom  is  assumed  independent  of  time.  When  these  assumptions 
are  made,  the  continuity  equation  (65)  takes  the  familiar  form: 

. |jOQ  = ° . (66) 
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2 . The  Keulegan  Lumped  Parameter  Model. 


The  basic  concept  of  the  limped  parameter  model  introduced  by  Keulegan 
(1967)  may  be  obtained  by  neglecting  the  right  side  of  equation  (64), 
assuming  that  ka  is  independent  of  x and  integrating  the  resulting 
equation  from  the  open  sea,  x = Xg,  to  the  bay,  x = Xg,  and  rearrang- 
ing terms  to  obtain: 


Z 


dy  dx  = 


f*  as 

J 3t 


(67) 


Keulegan  made  the  assumption,  consistent  with  equations  (4)  and  (5), 
that , 


(rzx^ Z 


K f U | U 

B2 


(68) 


and  implicitly  assumed  that  the  right  side  of  equation  (67)  could  be 
neglected. 


Substitution  from  equation  (68)  into  equation  (67)  yields: 


h - Ee]  ■ (Dpi  - 51]  * / 

l 

xs 

(69) 

~v.. 


Keulegan  (1967)  evaluated  uj,  not  at  the  same  value  of  x,  as  used  for 
the  tide  height  fig,  but  so  far  at  sea  that  a zero  value  could  be  assumed. 
He  implicitly  assumed  that  the  flow  was  uniform  throughout  the  channel. 
Thus,  the  velocity,  U/D,  is  given  by  ug.  Thus,  Keulegan's  basic  equa- 
tion takes  the  form: 


gfs  " ^g]  = ?*B 


yK(xB  ~ xsKo 

= un  • 

D B 


(70) 


The  two  terms  involving  the  square  of  the  speed  were  combined  to  give; 


1 - V 

2K(Xg  - xg) 

Lhs  hi 

d ‘ Y 

_ + m 

5 

_B 

2g 


(71) 


where  m is  a coefficient  related  to  the  velocity  profile  and  is  taken 
as  unity  in  most  applications j y = +1  for  floodflow,  and  y = -1  for 
ebb  flow.  Note  that  y(fig  - fig)  = |fig  - figl>  and  Y2  = 1. 
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Equation  (71)  contains  two  deviations  from  a rigorous  derivation, 
both  resulting  from  the  integration  of  the  term  3/3x  Q2.  A rigorous 
derivation  of  equation  (67)  requires  that  the  term  representing  the 
advection  of  momentum,  3u2/3x^  be  integrated  over  the  same  domain  as 
the  pressure  gradient  term,  3fi/3x,  so  that  the  gain  in  kinetic  energy 
associated  with  the  head  loss  resulting  from  the  slope  of  the  water  sur- 
face is  realized.  By  evaluating  the  water  level  at  opposite  ends  of  the 
inlet  and  the  kinetic  energy  in  the  open  sea,  where  u = 0,  and  at  the 
bay  end  of  the  inlet,  Keulegan  overestimated  the  gain  in  kinetic  energy 
as  a result  of  the  head  loss  within  the  inlet.  The  first  term,  u|  in 
equation  (70),  represents  kinetic  energy,  a positive  scalar  without  re- 
gard to  the  direction  of  current  flow.  TTie  second  term  represents  the 
bottom  stress,  a vector  which  changes  sign  with  changes  in  the  direction 
of  flow.  When  these  are  combined,  as  in  equation  (71),  the  kinetic 
energy  is  incorrectly  assigned  a negative  value  during  ebb  flow.  The 
kinetic  energy  term  is  often  much  smaller  than  the  friction  term  and  has 
only  a trivial  affect  on  the  solution  when  2K(x2  - Xj)/D  >>  m;  i.e., 
when  the  length  of  the  inlet  is  several  hundred  times  the  effective  depth. 
This  condition  is  satisfied  by  most  real  inlets. 

Keulegan  assumed  Manning's  formulation  of  the  bottom  friction  law 
where  K a D-1/3.  Since  this  assumption  is  not  critical  to  the  basic 
discussion  of  this  section,  K is  left  in  the  general  form,  K = K(x,  y) , 
without  considering  the  reason  for  variability. 

To  eliminate  the  velocity  uB  from  equation  (71),  Keulegan  assumed 
that  water  enters  and  leaves  the  bay  only  through  the  inlet,  the  area  of 
the  bay  is  independent  of  time,  and  the  phase  of  the  tide  in  the  bay  is 
independent  of  position  to  justify  the  expressions: 


or 


UBAB 


3hB 

= kBay  ~ * 


3t 


(72) 


where  V gay  and  ^Bay  are  t*ie  vol1™6  and  surface  area  of  the  bay.  These 
assumptions  are  rarely  satisfied  exactly,  but  may  be  approximated  with 
useful  accuracy  in  many  inlets. 


Elimination  of  u between  equations  (71)  and  (72)  yields: 


3h„ 


> 


(73) 
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where 


ab 

1 

OH 

ABay 

2K(xfi  - x^)  + mD 

(74) 


Ke  is  Keulegan's  coefficient  of  repletion,  and  apart  from  a few 
changes  in  notation,  equations  (73)  and  (74)  correspond  to  equations_12 
and  13  in  Keulegan  (1967).  All  five  empirical  parameters  Ag,  Ag^,  D, 

K,  and  (xg  -_Xg)  are  lumped  into  the  single  parameter,  Kg.  If  sufficient 
records  of  fig  and  h#  are  available  and  the  theory  holds,  can  be 

defined  in  a least  squares  sense  by  standard  statistical  procedures.  If 
the  various  assumptions  and  approximations  used  in  the  derivation  of  equa- 
tion (74)  are  approximately  valid,  it  should  be  possible  to  obtain  a good 
approximation  of  Kg  from  measurements  obtained  from  a hydrographic  chart. 


For  many  practical  problems,  the  assumptions  leading  up  to  equations 
(73)  and  (74)  are  too  restrictive  to  justify  much  confidence  in  a direct 
application  of  the  Keulegan  lumped  parameter  model.  The  principal  need 
for  a lumped  parameter  approach  is  to  provide  a method  for  estimating 
tidal  flows  in  small  or  dredged  inlets  where  observations  of  tides  are 
unavailable. 


Huval  and  Wintergerst  (1977)  propose  a relaxation  of  Keulegan's  (1967) 
assumptions  by  heuristic  extensions  of  his  basic  equations. 

The  differences  between  the  equations  used  and  those  found  in  a rig- 
orous derivation  were  not  identified  and  eliminated;  new  errors  were  intro- 
duced. To  display  these  errors,  their  importance  should  be  estimated,  and 
to  find  a satisfactory  solution,  it  is  perhaps  necessary  to  eliminate  the 
assumption  of  a constant  cross  section,  independent  of  x,  and  retain  the 
closure  terms.  The  one-dimensional  transport  equation  is  more  convenient 
for  this  purpose  than  the  mean  flow  equation  (64).  The  transport  equation 
CIO  in  Appendix  C,  has  the  form: 

y2 

+ ^x^c®2  + + / ^Tzx^Z  dy  = Rc  > (^s) 

where  Ac  = Ac(x)  = cross  section  of  inlet. 

3.  A Rigorous  Derivation. 

The  integral  of  equation  (75)  may  be  stated  as: 

= / [R*  ■ dx  • (76) 
xs 
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Terns  neglected  by  Keulegan  (1967)  are  on  the  right  in  equation  (76). 

The  first  term  in  this  equation  may  be  readily  integrated  to  obtain: 

rB  3 

J ax^112  dx  = ' (acu2)z5  • (77) 


If  either  of  the  factors,  A Q or  3h/3x  of  the  second  term  in  equa- 
tion (76),  is  nearly  constant,  the  integral  of  this  term  may  be  approxi- 
mated as: 

XB  - 

/ gAc|f  dx  = g(Rs  - %)  \Q  , (78) 

XS 

where  is  a representative  cross-sectional  area  for  the  entire  inlet. 

Equation  (78)  is  the  definition  of  Ae.  The  slope  of  the  water  surface, 
3fi/3y,  may  be  nearly  constant  whenever  the  phase  change  in  the  inlet  is 
small;  i.e.,  when 

(XB  ' xs0  <<:  [g Dmin1V2  T > C79) 


where  Dmin  is  the  minimum  value  of  D in  the  inlet,  and  T is  the 
period  of  the  tide  wave.  The  variability  of  Ae  can  be  determined  from 
a hydrographic  chart  or  a hydrographic  survey. 


Equation  (4)  may  be  used  to  obtain  the  integral  of  the  stress  term 
in  the  form: 

XB  y2  x0  y2 

/ / dy  dx  = y50  j j Ka)20*y)  dy  dx  > (80) 

XS  yl  XS  yl 

where 

w(x,y)  = iM/21  = . (81) 

UB  UB 


The  stress  coefficient,  K,  may  be  a function  of  x,  y,  and  t. 

Substitution  from  equations  (77),  (78),  and  (80)  into  equation  (76) 
yields : 


^BQB 


■ As^i 


♦ g(h5  - h5)  Ac  + yd 


u2 

'B 


XB  X2 

/ / 


Kw2  dy  dx 


3 

3t  c 


»] 


dx  . 


(82) 
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The  continuity  equation  for  the  one-dimensional  equations  of  motion 
is  derived  in  Appendix  C,  equation  (C9)  as: 


3A  e 
3t 


3x 


y2 

/ 


yi 


(R  + E)  dy  . 


(83) 


Equation  (75)  may  be  readily  integrated  along  the  inlet  to  obtain: 
XB 

8Ac 


/ 

XS 


at 


dx  + 


xb  y 2 

A55g  - A555  = j j (R  + E)  dy 


dx  . 


(84) 
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The  first  term  in  equation  (84)  may  be  simplified  by  noting  that 

XB 


j A^j  dx  = Vj.  , 


(85) 
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where  V|  is  the  volume  of  the  inlet.  Thus,  the  integrated  continuity 
equation  for  the  inlet  may  be  expressed  as: 

xr  y2 

avt  B 


- AS%  - A5u b + j j (R  + E)  dy 

xs  yi 


dx  . 


(86) 


When  the  phase  change  within  the  inlet  is  small,  the  volume,  as  a 
function  of  height,  may  be  computed  from  hydrographic  and  topographic 
data;  thus,  if  h5(t)  is  a known  function,  and  the  phase  change  within 
the  inlet  is  small,  Vj.(t)  may  be  determined. 


Equation  (86)  may  be  solved  for  u^  in  the  form: 

= AS=  1 3VI 

US  ■ 57*8  Aj  8t-  • 

where 

avT  av'  y2 

at-  " at  / / (R  + E)  dy  dx  . 

xs  Yl 


(87) 


(88) 


The  term,  aVj/at,  represents  the  contribution  of  tidal  flow  to 
the  local  time  derivative  of  the  inlet  volume.  The  term,  aVj/3t,  is 
the  complete  local  time  derivative  of  the  inlet  volume.  The  integral 
of  the  net  precipitation,  (R  + E) , over  the  area  of  the  inlet  is 
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the  net  contribution  of  precipitation,  evaporation,  and  overbank  flow 
to  the  volume  change  of  the  inlet.  Most  of  the  time  this  term  is  much 
smaller  than  the  volume  change  of  the  inlet  due  to  the  tidal  range  and 
can  be  neglected  in  most  cases.  An  exception  may  occur  during  intense 
rainstorms,  over  inlets  with  small  tidal  ranges.  In  most  cases,  3Vj/3t 
is  an  adequate  estimate  of  the  total  time  derivative  of  the  inlet  volume. 


Substitution  from  equation  (87)  into  equation  (82)  and  a slight  re- 
arrangement of  terms  in  equation  (70)  yield: 


Ag(h5  - hB) 


= 0 . 


(89) 


Equation  (89)  is  the  complete  form  of  the  momentum  equation  approxi- 
mated by  equation  (70).  It  is  a quadratic  equation  for  uB,  and  the 
solution  may  be  given  in  the  standard  form: 


= - b ± y/b^-  - 4 ac 

UB  = 2a 


(90) 


where 


a = y 


/ / Ku)2  dy  dx  + 
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It  is  necessary  to  neglect  the  second  term  of  a,  b,  and  the  last 
two  terms  of  c to  obtain  Keulegan's  (1967)  equation.  The  second  term 
in  a vanishes  for  a uniform  channel  cross  section  as  assumed  by  Keulegan 
but  the  other  terms  are  always  real  even  though  they  may  be  small  relative 
to  the  uncertainty  in  the  friction  term.  Consideration  of  typical  values 
for  each  term  in  equation  (90)  is  useful  to  provide  perspective  for  neglec- 
ting minor  terms,  remembering  that  the  initial  estimate  for  Manning's  n 
may  easily  be  in  error  by  50  percent.  Consequently,  K is  uncertain  by 
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a factor  of  two.  The  parameters  in  Table  1 may  be  taken  as  typical  of 
many  inlets  in  North  America.  In  an  order  of  magnitude  analysis,  inte- 
grals may  be  replaced  by  the  product  of  their  ranges  and  typical  values 
of  the  integrands.  Differentials  may  be  replaced  by  average  or  maximum 
values . 

Table  1.  Typical  inlet  parameters. 


Parameters 

Symbols 

Values 

Depth 

D 

3 m 

Width 

b 

300  m 

Cross-sectional  area 

9 x 102  m2 

Length 

LI 

2,000  m 

Surface  area 

b Lj 

6 x 10 5 m2 

Tidal  range 

2h 

1 m 

Volume  of  tidal  change 

2hbLj 

6 x 105  m3 

Tidal  period 

T 

12.4  hr 

Average  rate  of  inlet  volume  change 

3Vj/ 3t 

26.67  m3/ s 

Peak  velocity 

1 1 max 

1 m/s 

Manning's  n 

n 

0.025 

K 

n2D-« 

U) 

4.33  x 10“4 
1.0 

Acceleration  of  gravity 

g 

ab/a^ 

R - Rb 

10  m/s2 

0.9 

0.1m 

Substitution  for  most  terms  in  equation  (90)  yields: 

a = y(6  x 105)(4.33  x 10-4)  t 9 x 10  2 x 10"  1 
= 260  y ± 90 

b = 9 x lo_1  x 26.67 
= 24 

c = 9 x 102+1(h5  - hfl)  + 1.1  x 10"3  x 7.  nl  x 102 

+ 2 x 103  x [Rc  - 8 x 10“4] 

= 9 x 102  + 7.9  x 10'1  - 1.6  + 2 x 103RC 

It  is  difficult  to  assign  meaningful  typical  values  to  Rc.  Gener- 
ally, this  term  is  small  in  relation  to  the  first  term  in  c,  but  Rc 
may  be  the  largest  single  term  in  equation  (90)  when  a strong  wind  is 
blowing  parallel  to  the  inlet  axis.  The  wave  setup  term,  S^,  can  be 
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large  when  large  waves  are  breaking  in  the  inlet,  and  negligibly  small 
at  other  times.  Other  terms  may  also  become  significant  in  special  cir- 
cumstances . 


Substitution  of  numerical  values  into  equation  (89)  shows  that  for 
the  assumed  conditions,  3Vj/3t  and  3/3x  (Acu)  can  be  neglected  with- 
out serious  error,  even  when  average  absolute  values  are  considered. 
Since  these  are  alternating  terms,  the  actual  net  value  of  the  errors 
rising  from  neglect  is  generally  small.  The  second  term  of  a vanishes 
for  inlets  of  uniform  cross  section  as  assumed  by  Keulegan  (1967).  When 
the  cross  section  of  the  inlet  is  variable,  the  second  term  can  be  sig- 
nificant, and  for  short  inlets,  the  second  term  may  be  dominant.  Inlet 
dimensions  vary  widely,  and  values  for  the  particular  inlet  being  con- 
sidered should  be  substituted  into  equations  (89)  and  (90)  to  determine 
which  terms  can  be  neglected  in  particular  cases. 


Thus,  with  good  approximation,  ug  may  be  written  as: 

1 112 


us  = Y 


Ag|  h. 


■a  - "Bl 


*b  yi 

J j Ko)2  dy  dx  + yAg  (1  - Ag/A s) 
xs  yi 


(91) 


The  first  form  of  equation  (72)  may  be  generalized  to  account  for 
inflow  to  the  bay  from  precipitation,  overbank  flow,  and  rivers  by: 


3VD 

... Ba y = a A 

3t  BRBay 


(92) 


where  q^  is  the  net  contribution  from  rivers,  overbank  flow,  and  pre- 
cipitation, etc.  This  yields  an  expression  for  3hg/3t  in  the  form: 


% _ ^B = Hn 

3t  A Bay  B k Bay 


(93) 


Equation  (91)  can  be  used  to  eliminate  Ug  from  equation  (93)  to 
obtain : 


or 


3hg  _ YAg 

3t  kBay 


Ag|hg  - hg| 


XB  X2 

j j Ku2  dy  dx  + yAg(l  - Ag/Ag) 

yi 
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3 in 
kBay 


(94) 


(95) 
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where 


Ke  = 


A Bay 


XB  y2 

j j Km2  dy  dx  + yAb(1  - Ab/Ab) 


x5  yi 


V2 


(96) 


When  q^M  is  neglected,  equation  (95)  is  identical  with  equation 
(73).  Thus  equation  (95)  is  a generalization  of  equation  (73)  which 
permits  water  to  reach  the  bay  without  going  through  the  inlet.  When 
oj  = unity,  and  both  K and  A are  a constant,  equation  (96)  reduces  to 
the  form  of  equation  (74)  where  m is  neglected.  Thus,  equations  (95) 
and  (96)  are  generalizations  of  equations  (73)  and  (74). 

Huval  and  Wintergerst  (1977)  propose  a heuristic  generalization  of 
Keulegan's  (1967)  model  for  application  to  inlets  of  irregular  shape,  as 
partitioning  the  inlet  into  several  parallel  channels,  and  the  channels 
into  cells,  with  the  goal  of  obtaining  nearly  uniform  values  of  depth, 
roughness,  velocity,  and  cross-sectional  area  in  each  cell.  They  recog- 
nize that  these  conditions  cannot,  in  general,  be  satisfied  exactly; 
however,  no  guidance  is  provided  to  determine  how  well  these  conditions 
can  or  must  be  satisfied  to  justify  the  use  of  the  proposed  model. 


A systematic  approach  to  the  Huval  and  Wintergerst  model  leads  to  an 
alternate  expression  of  comparable  simplicity  which  does  not  require  such 
stringent  assumptions.  Moreover,  the  alternate  expression  can  be  con- 
verted into  the  form  given  by  Huval  and  Wintergerst  only  if  all  their 
assumptions  are  completely  satisfied. 


A rigorous  derivation  could  begin  with  the  numerical  integration  of 
equation  (Cl)  in  Appendix  C.  However,  the  principal  concern  here  is  the 
extension  of  Keulegan's  basic  equation  (70).  Thus,  the  derivation  can 
begin  with  a generalization  of  equation  (70),  based  on  equation  (69). 


It  will  be  assumed  that  an  inlet  can  be  partitioned  into  a number  of 
subchannels  and  that  within  each  subchannel,  the  depth,  roughness,  and 
velocity  are  nearly  constant  for  each  value  of  x.  In  general,  these 
variables,  ‘including  the  width,  within  each  subchannel  may  be  arbitrary, 
and  all  four  may  vary  with  x.  The  first  step  in  this  derivation  is  to 
rewrite  the  last  term  in  equation  (69)  as: 


k|u|u 

D2 


dy  dx  = 


Y 


Ku2 


dy  dx  . 


(97) 
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In  equation  (95)  , y j and  y2  are  taken  as  the  limits  of  the  sub- 
channel. With  the  aid  of  equation  (97)  and,  in  conformity  with  Keulegan 
(1967) , the  assumption  that  0^  can  be  neglected,  an  alternative  form 
of  equation  (70)  is  obtained: 


g[h,S  - hB]  = y 


f 

J A 


\n 


dx 


(9  8) 


The  subscript  m is  used  to  designate  the  m'th  subchannel. 

The  inertia  term,  arising  from  8Q2/8x,  has  been  omitted  as  explained 
above.  The  integration  of  equation  (98)  across  the  inlet  is  accomplished 
by  adding  the  solutions  for  all  subchannels  to  obtain: 

3 f 2 / ^P1  - h„)  . (99) 

m=  l fyn 

XS 


Equation  (81)  has  been  used  to  factor  5^  from  the  integral.  If  uB 
is  taken  as  the  average  velocity  across  the  entire  inlet,  equation  (99) 
may  be  solved  to  obtain: 


where 


_ 

Mg|R5  - hB| 

UB  = Y 

2 / K m<4/^  dx 

m=l  J 

L xs 

(100) 


M 

2 bm  = b , (101) 

m=i 


and,  as  in  equation  (91),  the  inertial  term  has  been  neglected.  Substi- 
tution from  equation  (100)  into  equation  (93)  yields  equation  (95),  with 
Ke  given  by: 


Ke  = 


B 

M 

&Bay 

M 

x5 

I 

f K dx 

m=i 

J 

XS 

— 

(102) 
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In  equations  (96)  and  (102),  the  major  problem  is  to  determine  the 
values  of  um(x) . Practical  procedures  are  developed  and  several  appli- 
cations of  the  procedures  are  given  by  Seelig,  Harris,  and  Herchenroder 
(in  preparation,  1977).  The  other  parameters  can  generally  be  measured 
or  estimated  from  available  charts. 


The  solution  of  Huval  and  Wintergerst  (1977)  may  be  obtained  from 
equation  (99)  by  setting  all  values  of 


Di 

so  that  equation  (99)  may  be  rewritten  as: 


(103) 


Mu 


5B  / 


LB  2 

Kl“? 


Di 


dx  = YgM(h5  - hs) 


(104) 


If  Ku2/D  is  the  same  function  of  x in  each  channel,  equation 
(98)  may  be  applied  separately  to  each  channel.  The  model  of  Huval  and 
Wintergerst  is  based  on  the  implicit  assumption  that  equation  (98)  can 
be  applied  separately  to  each  of  the  parallel  channels. 


Huval  and  Wintergerst  (1977)  defined  the  quantity,  w,  by  postulat- 
ing that  the  total  flow  in  each  subchannel  would  be  constant  for  all  x. 
By  using  this  assumption,  the  total  flow,  Q,  is  given  by: 

M 

Q = 2 % , (105) 

m=i 

where 

Qff?(x)  = %?(x)  » 

and 

M 

AB  = 2 » 

m=l 


substitution  from  equation  (81)  shows  that 


wm(x) 


i^(x)  - 

CB 


(106) 
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Factoring  equation  (104)  by  M,  substituting  from  equation  (106)  for 
uim,  and  rearranging  terms  in  equation  (99)  will  obtain: 


gl%  ~ hgl 

J dx/D^A^fx) 

xs 


(107) 


Summing  all  to  obtain  the  total  transport,  noting  that  Agug  = Q, 

and  substituting  the  sum  of  all  Qw  into  equation  (72)  yield: 


KHW 


1 M 

£ 

KBay  m=l 


I dx/DA^ (x) 


1/2 


(108) 


where  the  subscript,  HW,  is  used  to  indicate  the  form  developed  by  Huval 
and  Wintergerst. 

The  conditions  required  to  justify  equation  (108)  show  that  this  equa- 
tion is  equal  to  equation  (102).  All  terms  in  the  denominator  of  the 
bracketed  terms  in  equation  (108)  are  independent  of  y and  the  total 
cross  section  is  the  sum  of  the  cross  sections  of  the  individual  channels. 
If  A:  is  also  independent  of  x,  the  summation  in  equation  (108)  yields 


Kfftf  = 


A3 

1 

A Bay 

XB 

f Kj  dx/D 

4 

(109) 


For  the  conditions  required  to  justify  equation  (108),  i.e„,  Ax  is  con- 
stant and  K1  and  D are  independent  of  y,  all  terms  in  the  denomina- 
tor of  the  bracketed  term  in  equation  (102)  are  also  constant,  and  the 
total  width  of  the  inlet  is  equal  to  the  sum  of  the  widths  of  the  chan- 
nels. Thus,  for  these  conditions,  equation  (102)  may  be  written  as: 


where 


A B 

1 

kBay 

rB 

J dx/D 

_xs 

(HO) 


Aifl 

Al 


(x)  = 


= 1 


82 


The  general  validity  of  the  lumped  parameter  model  of  Keulegan  (1967) 
has  been  established  for  the  following  conditions: 

(a)  All  flow  into  the  bay  passes  through  the  inlet. 

(b)  Surface  area  of  inlet  and  bay  is  independent  of  time. 

(c)  Cross-sectional  area  of  inlet  is  independent  of  x. 

(d)  Phase  change  of  the  tide  within  inlet  and  bay  is  small. 

(e)  The  friction  effect  is  guch  larger  than  the  effect  of  momentum 
advection;  e.g.,  2K(xs  - xg)  >>  5. 

This  latter  condition  is  required  to  ensure  that  the  error  in  treat- 
ing the  terms  arising  from  the  term  3u2/8x  does  not  become  significant. 

Huval  and  Wintergerst  (1977)  successfully  relax  the  following  condi- 
tions : 

(a)  Full  relaxation. 

(b)  Area  of  inlet  and  bay  may  change  as  linear  or  quadratic  func- 
tions of  R.  This  is  sufficiently  generalized  for  most  practical  cases. 

(c)  The  requirement  for  a channel  cross  section  independent  of  x 
is  generalized  as  follows: 

(1)  The  inlet  is  partitioned  into  several  channels  of  equal 
width.  The  width  may  change  with  distance  along  the  inlet  axis,  but  the 
cross  section  in  each  channel  is  as  nearly  independent  of  x as  possible. 

(2)  The  transport  in  each  channel  is  constant  along  its  length. 
This  is  equivalent  to  the  specification  that  the  subchannel  boundaries  be 
streamlines . 

(3)  The  factor,  has  the  same  value  in  each  channel 

for  each  value  of  x. 

Condition  (c)  is  unlikely  to  be  satisfied  even  approximately  unless 
the  inlet  is  approximately  symmetric  to  the  centerline.  Conditions  (d) 
and  (e)  for  the  Keulegan  (1967)  model  are  retained  in  the  model  of  Huval 
and  Wintergerst  (1977). 

Both  models  become  invalid  whenever  strong  winds  blow  parallel  to  the 
inlet  axes  and  when  large  waves  break  in  the  inlet  There  may  be  special 
but  unusual  circumstances  when  either  of  the  models  fails  because  one  or 
more  other  processes  grouped  under  the  term  Ra  become  important. 
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4.  Summary. 


Equations  of  motion  for  the  vertically  and  cross-channel  integrated 
flow  in  a channel  of  variable  cross  section  have  been  derived  in  Appendix 
C in  a manner  which  explicitly  identifies  the  assumptions  generally  made 
in  adopting  one-dimensional  equations  of  motion.  These  equations  have 
been  further  integrated  along  the  channel  in  this  section  to  provide  a 
differential  equation  for  the  tidal  height  in  the  ocean  and  the  inlet 
geometry.  The  work  was  motivated  by  the  inlet  studies  of  Keulegan  (1967) 
and  Huval  and  Wintergerst  (1977). 

This  study  has  shown  that  the  inertia  terms  appear  to  have  been  im- 
properly treated  by  Keulegan  (1967) ; however,  these  terms  are  of  trivial 
importance  in  some  inlets  and  may  have  little  affect  on  the  flow. 

Huval  and  Wintergerst  attempt  to  generalize  Keulegan's  treatment  by 
relaxing  some  of  his  restriction?,  without  going  back  to  a complete  set 
of  equations  as  a starting  point.  A thorough  analysis  of  the  problem 
has  shown  that  Huval  and  Wintergerst  found  a singular  solution  of  the 
problem  which  is  nearly  as  restrictive  as  that  presented  by  Keulegan. 

The  new  derivation  in  this  study  produced  a solution  of  the  same  level 
of  complexity  which  agrees  with  Huval  and  Wintergerst  whenever  their 
solution  is  valid,  and  which  is  valid  for  a much  wider  range  of  condi- 
tions. The  new  solution  has  not  been  developed  here  for  practical  appli- 
cation, but  is  accomplished  in  Seelig,  Harris,  and  Herchenroder  (in  pre- 
paration, 1977). 

VI.  MODEL  CALIBRATION 

Both  the  numerical  and  physical  models  used  for  the  study  of  coastal 
hydrodynamics  are  designed  to  reproduce,  as  reliably  as  possible,  the  pri- 
mary flow  phenomena  being  investigated,  and  to  approximate  or  simulate  in- 
directly the  effects  of  secondary  phenomena.  The  optimum  design  simula- 
tion and  approximation  procedures  cannot  be  determined  or  evaluated  solely 
on  the  basis  of  theory.  Thus,  a process  of  experimental  adjustment  and 
evaluation  of  the  model  must  be  carried  out  before  it  can  be  used  with 
confidence  for  either  engineering  or  scientific  studies.  This  process 
has  been  termed  confirmation , adjustment , verification , or  calibration 
by  various  investigators.  In  this  study,  calibration  is  used  to  describe 
the  process  where  a physical  or  numerical  model  of  hydraulics  phenomena 
is  checked  with  prototype  data  and  systematically  adjusted  to  reproduce 
desired  prototype  data  from  corresponding  prototype  driving  forces.  Ver- 
ification is  restricted  to  the  process  where  independent  prototype  data 
(data  not  used  in  the  calibration)  are  used  to  verify  that  a calibrated 
hydraulics  model  reproduces  prototype  observations  with  satisfactory 
accuracy. 

Models  of  coastal  phenomena  can  be  useful  without  being  completely 
accurate.  Hudson,  et  al.  (in  preparation,  1977,  pp.  3-5)  stated  that 
all  models  of  estuaries  have  one  characteristic  in  common:  that  is  that 
they  cannot  be  completely  accurate  simulations  of  all  the  complex  phe- 
nomena inherent  in  tidal  waterways....” 
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American  Society  of  Civil  Engineers  (1942)  reported  "That  the  neces- 
sity for  the  trial- and-error  process  of  model  verification  arises  from 
the  inability  of  the  experimenter  to  design  and  construct  his  model  so 
as  to  insure  the  attainment  of  dynamic  similitude."  The  previous  quota- 
tion (Hudson,  et  al.,  in  preparation,  1977)  showed  that  this  need  is  still 
much  in  evidence. 

American  Society  of  Civil  Engineers  (1942)  and  Hudson,  et  al.  (in 
preparation,  1977)  give  lengthy  descriptions  of  the  procedures  used  for 
the  adjustment  of  physical  models.  Durham,  Greer,  and  Whalin  (1976)  and 
Whalin,  Perry,  and  Durham  (in  preparation,  1977)  discuss  a calibration 
procedure  based  on  the  harmonic  analysis  of  tidal  records  for  a period 
of  about  15  days  which  appears  to  be  an  improvement  on  traditional  pro- 
cedures. The  field  data  collection  program  at  Masonboro  Inlet  was  too 
restricted  to  support  this  approach.  Odd  (1971)  surveyed  the  problem 
in  less  detail  for  one-dimensional  mathematical  models.  The  four  appen- 
dixes to  this  report  present  specific  examples  of  model  calibration.  This 
section  does  not  furnish  a guide  to  model  calibration,  but  provides  the 
reader  with  an  understanding  of  the  physical  and  mathematical  need  for 
calibration  and  with  background  information  needed  for  the  assessment  of 
the  calibration  procedures  for  the  four  Masonboro  Inlet  models. 

The  imperfections  in  hydraulics  models  that  lead  to  the  requirement 
for  calibration  result  in  part  from  the  lack  of  precise  information  about 
the  prototype  and  in  part  from  distortion  or  total  neglect  in  the  model 
of  processes  that  are  of  secondary  but  significant  importance  in  the  pro- 
totype. 

Friction  is  generally  the  most  important  of  the  distorted  physical 
processes.  The  process  of  model  adjustment  for  improving  the  calibration 
of  model  to  prototype  is  often  discussed  as  if  the  only  important  problem 
is  to  secure  the  proper  value  of  bottom  friction  everywhere  in  the  model. 

Friction  has  two  important  affects  on  the  prototype  flow.  Turbulence., 
generated  by  friction,  serves  to  mix  adjacent  fluid  elements,  thus  smooth- 
ing out  large-scale  gradients  of  momentum  and  dispersing  any  substance 
mixed  with  the  water.  Turbulence  also  increases  the  frequency  of  locally 
intense  gradients  in  the  flow,  thereby  increasing  the  conversion  of  kinetic 
energy  to  heat  by  molecular  viscosity.  These  two  effects  of  turbulence 
and  one  other  effect  are  important  in  the  physical  model.  By  increasing 
the  generation  of  turbulence  near  the  bottom,  friction  decreases  the 
thickness  of  the  viscous  boundary  layer,  and  the  thickness  of  the  fluid 
layer  in  which  molecular  viscosity  exerts  a direct  control  over  the  veloc- 
ity profile  in  the  fluid. 

The  bottom  stresses  must  be  increased  by  a factor  equal  to  the  distor- 
tion ratio,  B/a,  to  be  correctly  represented  in  distorted  scale  physical 
models  (see  Sec.  III).  Hudson,  et  al.  (in  preparation,  1977)  derived  the 
same  correction  factor  by  a different  procedure. 

Internal  friction,  as  described  by  all  nine  components  of  the  turbu- 
lent stress  tensor  (App.  B,  eq.  B13),  can  contribute  to  fluid  mixing  and 
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the  dissipation  of  energy.  When  no  special  effort  is  made  to  represent 
any  one  of  these  terms  in  the  model,  whether  physical  or  numerical,  the 
missing  processes  are  assumed  to  be  lumped  with  the  bottom  stresses  in 
the  model.  Although  not  always  clearly  expressed,  this  principle  is 
often  used  in  the  adjustment  of  physical  models.  The  most  frequently 
used  procedure  at  the  U.S.  Army  Engineer  Waterways  Experiment  Station 
(WES),  Vicksburg,  Mississippi,  is  to  embed  small  metal  strips  into  the 
concrete  bottom  of  fixed-bed  models  and  to  bend  or  remove  the  strips 
later  as  required  to  obtain  a satisfactory  level  of  turbulence  during 
the  model  calibration  procedure.  This  procedure  and  others  are  described 
by  Sager  and  Seabergh  (19  77). 

However,  according  to  Hudson,  et  al.,  (in  preparation,  1977),  it  is 
not  always  possible  to  achieve  sufficient  mixing  of  fluid  elements  by 
these  techniques,  and  oscillating  fans  may  be  directed  toward  the  water 
surface  to  increase  the  intensity  of  turbulence  in  the  water.  Small 
electrical  vortex  generators  similar  to  kitchen  mixers  are  used  in  some 
European  laboratories  for  the  same  purpose. 

The  geometry  of  the  inland  tide  reservoirs  (see  Sec.  Ill),  may  be 
changed  to  improve  the  calibration  of  the  physical  model.  Sometimes  it 
is  necessary  to  deepen  a noncritical  part  of  the  basin  or  to  constrict 
the  flow  somewhere  in  the  model  more  than  strict  geometric  similitude 
would  prescribe,  to  duplicate  some  important  features  of  the  prototype 
flow  in  the  model. 

Reproduction  of  water  levels  is  more  easily  achieved  than  the  repro- 
duction of  currents  in  both  physical  and  mathematical  models.  Hudson, 
et  al.  (in  preparation,  1977)  recommended  that  acceptable  calibration  of 
water  levels  through  one  or  more  tidal  cycles  be  obtained  before  attempt- 
ing to  adjust  the  model  for  the  reproduction  of  currents. 

Model  adjustments  which  are  not  determined  by  theoretical  considera- 
tions leave  something  to  be  desired,  but  are  often  considered  necessary. 

A model  that  cannot  reproduce  recorded  events  acceptably  cannot  be  de- 
pended on  to  predict  future  events.  The  final  stage  of  calibration  for 
many  models  appears  to  be  a trial -and- error  procedure. 

Numerical  hydraulic  models  have  not  been  in  use  as  long  as  physical 
models.  Consequently,  the  techniques  for  adjustment  during  the  calibra- 
tion process  are  not  well  developed.  In  general,  the  development  has 
followed  along  the  lines  previously  used  with  physical  models.  Some 
techniques,  however,  are  more  conveniently  applied  to  one  type  of  model 
or  the  other. 

The  concept  of  considering  all  adjustments  of  the  model  as  changes  in 
the  bottom  stress  has  carried  over  from  physical  modeling.  In  numerical 
models,  the  stress  can  easily  be  expressed  in  the  form  of  equations  (4) 
or  (5)  with  a constant  stress  coefficient.  Also,  an  initial  estimated 
stress  coefficient  which  may  vary  with  position,  can  easily  be  changed 
by  a constant  factor  or  a constant  increment.  Since  this  type  of  stress 
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law  or  change  in  the  stress  law  is  difficult  to  achieve  in  a physical 
model,  it  is  used  more  often  with  numerical  models.  Reid  and  Bodine  (1968) 
used  a constant  stress  coefficient  for  the  entire  model  of  Galveston  Bay, 
but  varied  the  value  between  repeated  calculations  to  obtain  optimum 
agreement  between  model  and  prototype.  Chiu,  van  de  Kreeke,  and  Dean 
(1970)  also  experimented  with  constant  values  for  K during  each  cali- 
bration and  selected  the  value  which  optimized  agreement  between  model 
and  prototype  for  later  experiments.  Masch,  Brandes,  and  Reagan  (1977) 
chose  initial  values  of  K which  varied  with  water  depth,  and  later 
adjusted  some  values  on  a point  by  point  basis  to  optimize  the  agreement. 
Chen  and  Hembree  (1977)  intended  to  follow  a similar  procedure;  their 
computer  program  has  a provision  for  reading  in  an  initial  estimate  of 
friction  factors.  However,  the  input  friction  factors  are  not  included 
with  the  program  listing,  and  no  numerical  values  or  descriptive  logic 
are  provided.  Thus,  it  is  impossible  to  determine  from  the  material 
submitted,  how  the  bottom  stress  was  actually  treated  in  this  numerical 
model . 

Odd  (1971)  reviewed  the  problem  of  representing  bottom  friction  in 
one-dimensional  numerical  tidal  models  and  concluded  that  empirical  re- 
sistance functions,  such  as  the  Chezy  and  Manning  formulas,  provided  an 
inadequate  description  of  frictional  resistance  in  shallow  estuaries  with 
an  appreciable  tidal  range.  The  reasons  for  this  opinion  are  not  pre- 
sented with  enough  details  for  application  to  two-dimensional  models. 

Reid  and  Bodine  (1968)  treated  the  overtopping  of  previously  dry  grid 
points  by  imposing  a simulated  weir  and  supplying  both  weir  elevations 
and  coefficients  for  each  case.  Thus,  additional  constants  become  avail- 
able for  adjustment  during  the  calibration  process.  Simulated  weirs  could 
be  used  to  impose  a reduction  in  the  cross  section  of  flow  to  compensate 
for  the  difficulty  in  representing  channels  with  widths  less  than  As, 
or  channels  which  could  not  otherwise  be  properly  represented  within  the 
selected  network  of  grid  points.  Procedures  similar  to  those  used  by 
Reid  and  Bodine  (1968)  were  applied  in  both  of  the  Masonboro  Inlet  two- 
dimensional  numerical  models. 

None  of  the  reports  cited  describe  the  processes  used  in  deciding  how 
roughness  elements  should  be  distributed  in  a physical  model,  or  how  quan- 
titative coefficients  should  be  selected  in  numerical  models  to  efficiently 
improve  the  agreement  between  model  and  prototype.  Therefore,  model  cali- 
bration by  both  physical  and  numerical  models  remains  a trial-and-error 
process . 

Adjustment  of  a model  to  optimize  agreement  between  the  model  and 
prototype  for  known  events,  and  verification  to  measure  the  degree  of 
success  achieved,  are  necessary  steps.  A model  cannot  be  depended  on  to 
reproduce  future  events  with  greater  accuracy  than  the  model  has  displayed 
in  reproducing  known  events.  However,  there  is  no  guarantee  that  the  pre- 
diction of  future  events  will  achieve  the  degree  of  success  displayed  in 
calibration  even  when  the  model  is  not  modified.  Any  substantial  modifi- 
cation of  the  model  is  likely  to  change  the  calibration. 
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The  reliability  of  both  physical  and  numerical  models  for  predicting 
the  effects  of  changed  conditions  will  be  greatly  increased  when  the 
process  of  adjustment  is  changed  from  an  art  to  a science. 

VII.  COMPARISON  OF  MODEL  RESULTS 

The  performance  of  each  model  used  in  the  Masonboro  Inlet  study  for 
simulating  tidal-induced  water  motions  is  discussed  in  this  section.  The 
inlet  modeling  study  consisted  of  two  phases:  (a)  Calibration  of  the 
models  to  1969  basin  geometry  and  to  given  1969  velocity  and  tidal  data, 
and  (b)  application  of  the  models  to  Masonboro  Inlet  for  basin  geometry, 
corresponding  to  prejetty  (1964)  conditions.  In  the  first  phase,  the 
effectiveness  of  each  model  in  reproducing  water  motions  is  shown  by  com- 
paring the  tides  and  velocities  predicted  in  the  models  to  the  tides  and 
velocities  observed  in  the  1969  prototype;  in  the  second  phase,  the  re- 
sults obtained  with  the  mathematical  models  are  compared  with  the  tides 
and  velocities  measured  in  the  physical  model  for  1964  conditions. 

1.  Calibration  of  Models. 


Basic  data  for  calibrating  the  models  to  Masonboro  Inlet  were  from 
basin  geometry  collected  during  September  and  October  1969,  and  from 
tidal  excitations  and  current  velocities  collected  on  12  September  1969. 
Basin  geometry  data,  including  both  hydrography  and  topography,  were 
taken  from  a chart  (scale  of  1 inch  to  500  feet)  with  2-foot  (0.61  meter) 
contours  representing  the  hydrography  below  the  mean  low  water  (MLW) 
Beaufort  Datum,  or  1.88  feet  (0.51  meter)  below  the  National  Geodetic 
Vertical  Datum  (formally  called  the  1929  Mean  Sea  Level  Datum)  elevations 
for  describing  the  ground  elevations  above  the  local  mean  high  water  (MHW) 
Beaufort  Datum.  Selection  of  the  physical  limits  of  water  motions  at 
specific  locations  was  reproduced. 

Tide  and  current  velocity  data  were  collected  at  several  locations. 

For  the  period  0300  to  2100  hours,  eastern  standard  time  (e.s.t.),  con- 
tinuous tidal  information  was  recorded  by  seven  tide  gages.  Concurrently, 
from  about  0600  to  2100  hours,  e.s.t.,  current  velocities  were  measured 
and  recorded  at  15  stations.  Locations  of  the  tide  gages  and  velocity 
measurement  stations  are  shown  in  Figure  1;  three  velocity  measurement 
stations  were  used  to  measure  the  velocities  in  a particular  channel 
cross  section.  Measurements  were  taken  at  each  station  at  three  levels 
in  the  vertical,  1 foot  above  the  bottom,  middepth,  and  1 foot  below  the 
free  surface. 

Tides  recorded  on  12  September  1969  at  gage  0,  located  on  Mercer's 
Pier  on  the  Oceanside  of  Wrightsville  Beach  and  about  15  miles  from 
Masonboro  Inlet  (Fig.  1),  provided  the  basic  input  to  all  models.  The 
data  recorded  at  gage  0 (Fig.  4)  indicates  that  the  tide  for  this  loca- 
tion is  nearly  sinusoidal  in  shape  and  has  a semidiurnal  period  of  about 
12.5  hours. 
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2.  Comparison  of  Results. 


Tides  reproduced  in  the  two-dimensional  mathematical  models  and  physi- 
cal model  in  comparison  to  the  tides  recorded  on  12  September  1969  at  the 
various  gage  locations  are  shown  in  Figures  5 to  10.  All  plotted  tidal 
elevations  were  referred  to  the  MLW  Beaufort  Datum.  Tides  reproduced  in 
both  the  physical  model  (Sager  and  Seabergh,  1977)  and  in  the  numerical 
model  by  Masch,  Brandes,  and  Reagan  (1977)  used  a datum  level  that  was 
intended  to  correspond  with  the  MSL.  However,  each  model  used  a different 
constant  for  converting  from  the  MLW  datum  to  the  MSL  datum.  In  each  case, 
the  respective  conversion  factor  was  used  to  convert  the  tidal  elevations 
back  to  the  MLW  datum. 

The  ocean  tide  recorded  at  gage  0 is  superposed  on  the  tide  hydro- 
graphs in  Figures  5 to  10  to  demonstrate  the  modification  of  the  tidal 
curves  with  respect  to  both  phase  and  waveform  as  the  tide  propagates 
through  the  entrance  channel  and  into  the  back-bay  channels.  Only  a 
slight  change  occurs  in  phase  and  waveform  in  the  entrance  channel  (Fig. 

5).  As  expected,  the  wave  in  the  inlet  and  back-bay  channel  lags  more  in 
phase,  and  the  waveform  becomes  more  distorted.  Generally,  the  tide  wave 
is  not  greatly  modified  by  the  physical  characteristics  of  the  Masonboro 
Inlet  study  area.  The  root-mean-square  (rms)  difference  between  proto- 
type and  model  water  levels  varied  from  0.1  to  0.3  foot  for  individual 
gages,  with  an  overall  value  of  0.2  foot  for  each  model.  None  of  the 
models  appeared  significantly  superior  to  the  others  for  reproducing 
the  tides  at  all  tidal  stations. 

Figure  11  shows  the  tide  predicted  for  the  embayment  area  with  the 
lumped  parameter  model  in  comparison  to  the  average  embayment  tide 
observed  in  the  prototype.  The  average  embayment  tide  for  the  prototype 
was  calculated  by  averaging  the  tides  observed  at  gage  stations  3,  4, 
and  5.  This  average  was  used  because  the  lumped  parameter  mathematical 
model  can  only  predict  a uniform  variation  of  the  water  surface  in  the 
embayment  area.  Although  the  figure  indicates  that  the  amplitude  of  the 
tide  calculated  at  the  maximum  stage  was  too  large,  the  predicted  tide  is 
generally  in  good  agreement  with  the  observed  tide. 

A comparison  of  the  predicted  and  prototype  current  velocities  for 
the  calibration  data  of  12  September  1969  is  shown  in  Figures  12  to  26. 
Surface  measurements  were  made  0.30  meter  (1  foot)  below  the  surface; 
bottom  measurements  were  made  0.30  meter  above  the  bottom.  The  proto- 
type velocities  are  based  on  measurements  taken  from  boats  using  Price 
cuptype  velocity  meters.  The  instrumental  error  is  less  than  10  percent 
when  the  recorded  speed  exceeds  3 centimeters  (0.1  foot)  per  second. 

Velocities  measured  in  the  physical  model  were  taken  in  the  vertical 
at  each  station  similar  to  the  prototype  (Figs.  12  to  26).  The  terms 
surface  and  bottom  in  the  figures  are  interpreted  as  1 foot  below  the 
water  surface  and  1 foot  above  the  bed,  respectively.  Unlike  the  physi- 
cal model,  the  velocities  calculated  in  the  mathematical  models  are 
vertically  averaged  velocities. 
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Figure  5.  Observed  and  predicted  tides  at  gage  1,  12  September  1969. 
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Figure  7.  Observed  and  predicted  tides  at  gage  3,  12  September  1969. 


LEGEND 


Observed  and  predicted  tides  at  gage  5,  12  September  1969. 
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Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  1,  station  C (center) 
12  September  1969. 
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Figure  14.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  1,  station  S (south) 
12  September  1969. 
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Figure  19.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  3,  station  C (center) , 
12  September  1969. 
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Figure  20.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  3,  station  S (south), 
12  September  1969. 
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Figure  22.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  4,  station  C (center), 
12  September  1969. 
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Figure  23.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  4,  station  E (east), 
12  September  1969. 
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Figure  26.  Prototype  and  predicted  velocity  for  Masonboro  Inlet  at  range  5,  station  W (west), 
12  September  1969. 


The  velocity  curves  (Figs.  12  to  26)  show  that  Chen  and  Hembree 
(1977)  were  unable  to  achieve  a good  velocity  verification.  Although 
the  phasing  was  in  fairly  good  agreement  with  observations,  the  magni- 
tude of  the  flood  and  ebb  velocities  was  generally  poor  with  respect  to 
the  prototype  velocities.  However,  both  Sager  and  Seabergh  (1977)  and 
Masch,  Brandes,  and  Reagan  (1977)  obtained  a better  velocity  verification. 

The  lumped  parameter  mathematical  model,  because  of  its  formulation, 
can  predict  only  the  mean  velocities  in  the  throat  of  the  inlet.  To  show 
how  well  this  model  was  calibrated  to  prototype  conditions,  reduction  of 
the  velocities  measured  in  the  inlet  to  mean  velocities  was  necessary. 

To  obtain  mean  velocities  for  the  inlet  from  the  measurements,  all  nine 
measurements  for  the  three  velocity  stations  were  averaged.  A simple 
average  was  used  to  determine  both  the  predicted  and  prototype  mean  ve- 
locities. A comparison  of  the  inlet  velocities  reproduced  in  the  lumped 
parameter  model  and  the  mean  prototype  velocities  is  shown  in  Figure  27. 
The  figure  also  shows  the  mean  velocities  determined  by  Masch,  Brandes, 
and  Reagan  (1977)  and  Sager  and  Seabergh  (1977).  Results  are  not  shown 
in  Chen  and  Hembree  (1977)  because  of  the  poor  velocity  reproductions. 

Both  the  physical  and  numerical  models  by  Masch,  Brandes,  and  Reagan 
simulated  the  prototype  velocities  at  all  stations  more  satisfactorily 
than  the  model  by  Chen  and  Hembree  (see  Figs.  12  to  26).  However,  it  is 
difficult  to  determine  whether  the  physical  model  or  the  numerical  model 
by  Masch,  Brandes,  and  Reagan  gave  the  best  overall  reproduction.  A 
systematic  procedure  for  making  this  determination  is  to  compare  the 
rras  difference  between  prototype  velocities  and  the  velocities  obtained 
from  the  models. 

To  carry  out  such  an  analysis,  the  prototype  and  simulated  velocities 
were  tabulated  hourly  for  each  station.  Velocities  measured  in  the  phys- 
ical model  and  prototype  were  converted  to  vertically  averaged  velocities 
to  agree  with  the  velocities  reproduced  in  the  numerical  model.  These 
conversions  were  made  by  taking  a simple  average  of  the  velocities,  usu- 
ally given  at  three  levels  in  the  vertical,  for  a given  time  level.  The 
root-mean-square  velocity  difference  (RMSA)  was  determined  for  each 
station  in  both  models.  The  RMSA  is  defined  by: 

RMSA  = [i  Z(Y  - X)2f  , (Hi) 

where  Y is  the  model  velocity,  X is  the  prototype  velocity  for  the 
same  observation  time,  and  N is  the  number  of  variables  compared.  Hie 
RMSA's  were  normalized  by  division  by  the  mean-square  value  of  the  pro- 
totype measurement  RMSP,  defined  by: 

RMSP  = EX2] V2  . (112) 
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The  results  of  this  analysis  for  each  velocity  station  are  shown  in 
Table  2.  The  last  two  columns  in  the  table  give  the  relative  normalized 
variations,  N( RMSA/RMSP)2,  for  the  physical  and  Masch,  Brandes,  and 
Reagan  (1977)  models.  A comparison  of  the  values  in  the  columns  implies 
that  velocities  are  reproduced  more  satisfactorily  at  some  stations  in 
the  physical  model  and  at  other  stations  in  the  model  of  Masch,  Brandes, 
and  Reagan.  If  each  of  the  last  two  columns  in  Table  2 is  summed  and 
divided  by  the  sum  of  N,  the  overall  normalized  mean  squares  are  found. 
The  square  root  of  the  normalized  mean  squares  is  simply  the  normalized 
root  mean  squares,  RMSA/RMSP.  By  carrying  out  these  calculations,  it  is 
found  that  RMSA/RMSP  = 0.355  for  the  physical  model  and  0.299  for  the 
model  by  Masch,  Brandes,  and  Reagan.  If  those  stations  are  omitted  when 
velocity  data  were  missing  at  some  levels  in  the  vertical  in  the  physical 
model,  then  RMSA/RMSP  = 0.313  in  the  physical  model  and  0.300  in  the  model 
by  Masch,  Brandes,  and  Reagan.  Since  only  slight  differences  in  the  over- 
all rms  deviations  were  found,  it  is  impossible  to  conclude  that  either 
model  reproduced  the  velocities  more  satisfactorily  than  the  other  because 
of  possible  errors  in  reading  the  velocities  from  curves  and  only  a lim- 
ited set  of  data  was  used  in  the  analysis.  In  general,  both  models  per- 
formed about  equally  well  in  reproducing  the  velocities.  However,  it 
should  not  be  assumed  that  the  best  possible  calibrator  was  obtained  for 
any  of  the  models. 

3.  Application  of  the  Models. 

The  second  phase  of  the  modeling  study  was  applying  the  calibrated 
models  of  Masonboro  Inlet  to  calculations  for  periods  where  the  geometry 
of  the  system  was  substantially  different  from  that  for  the  1969  condi- 
tions. Preproject  conditions  (1964)  when  Masonboro  Inlet  was  in  a natural 
state  were  selected. 

Because  tides  or  velocities  were  not  available  for  the  preproject  con- 
ditions, idealized  mean  and  spring  tides  characteristic  of  the  Masonboro 
Inlet  area  were  used  for  excitation  of  the  models.  The  idealized  tidal 
curves  are  shown  in  Figure  28. 

Application  of  the  models  to  Masonboro  Inlet  for  the  preproject  con- 
ditions was  made  primarily  to  obtain  additional  insight  on  the  perform- 
ance of  numerical  models.  If  the  initial  calibration  of  the  physical 
model  had  been  significantly  better  than  either  numerical  model,  and  if 
the  divergence  between  physical  and  numerical  models  had  increased  sig- 
nificantly for  the  preproject  conditions,  the  assumption  would  be  that 
neither  of  the  numerical  models  considered  would  be  reliable  for  appli- 
cation to  any  condition  other  than  that  used  in  the  initial  calibration. 
Due  to  the  small  differences  between  the  idealized  mean  and  spring  tides, 
the  comparison  of  the  modeling  result  is  shown  only  for  the  mean  tide  con- 
ditions (Fig.  28).  Information  was  furnished  for  modifying  the  models  to 
reflect  the  hydrography  corresponding  to  the  November  1964  preproject  con- 
ditions, (Fig.  29). 
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Figure  28.  Mean  and  spring  tides  at  Masonboro  Inlet. 
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Figure  29.  Hydrographic  survey  for  Masonboro  Inlet,  November  1964. 


The  variations  in  water  surface  elevations  predicted  by  all  three 
two-dimensional  models  are  shown  by  hydrographs  in  Figures  30  to  35. 

These  hydrographs  show  that,  in  general,  the  tides  simulated  by  the 
two-dimensional  mathematical  models  were  in  reasonable  agreement  with 
those  simulated  in  the  physical  model.  However,  about  a %-hour  phase 
lag  occurs  in  the  tides  predicted  by  Masch,  Brandes,  and  Reagan  (1977) 
relative  to  the  predictions  by  the  physical  model  at  gage  6.  Chen  and 
Hembree  (1977)  did  not  calculate  the  tide  for  this  location.  Since  the 
tides  calculated  by  Masch,  Brandes,  and  Reagan  matched  the  prototype 
observations  at  station  6 better  than  the  tides  observed  in  the  physical 
model,  the  numerical  model  possibly  provided  a better  approximation  than 
the  physical  models  on  the  tides  at  this  station. 

Figure  36  shows  the  uniform  embayment  water  level  elevations  calcu- 
lated by  the  lumped  parameter  model  in  comparison  to  the  average  embay- 
ment elevations  reproduced  in  the  physical  model;  the  lumped  parameter 
model  is  in  good  agreement  with  the  physical  model  in  reproduction  of 
tides  for  the  preproject  conditions. 

Predicted  tidal  current  velocities  for  all  three  models  are  shown 
in  Figures  37  to  51.  The  agreement  between  the  physical  model  and  the 
numerical  model  by  Masch,  Brandes,  and  Reagan  is  again  much  better  than 
the  agreement  between  either  of  these  models  and  the  model  by  Chen  and 
Hembree.  However,  there  are  significant  differences  between  the  physical 
model  and  the  numerical  model  by  Masch,  Brandes,  and  Reagan.  Stronger 
ebb  current  velocities  were  computed  by  Masch,  Brandes,  and  Reagan  at 
station  3C  (Fig.  44);  smaller  ebb  current  velocities  were  computed  by 
this  model  at  stations  4E,  4W,  5C,  and  5W  (Figs.  46  and  48  to  51).  Why 
these  deviations  occurred  or  which  prediction  is  more  correct,  is  unknown. 

The  velocities  computed  by  the  lumped  parameter  method  (Huval  and 
Wintergerst,  1977)  are  compared  to  average  velocities  in  the  inlet  as 
predicted  by  both  the  physical  model  (Sager  and  Seabergh,  1977)  and  the 
numerical  model  (Masch,  Brandes,  and  Reagan,  1977)  (Fig.  52). 

4.  Later  Studies. 


A new  hydrographic  survey  of  Masonboro  Inlet  was  made  in  July  1974. 
The  physical  model  of  the  inlet  was  remodeled  to  the  1974  bathymetry  and 
the  measurements  described  above  were  repeated.  The  rms  differences  be- 
tween model  and  prototype  and  two  other  measures  of  goodness  of  calibra- 
tion were  computed.  The  results  are  given  by  Seabergh  and  Mason  (1975) . 
A report  covering  both  physical  and  numerical  models  will  be  published 
later  under  the  GITI  program.  The  rms  differences  were  slightly  larger 
than  in  calibration  tests,  as  would  be  expected  when  no  effort  was  made 
to  force  the  model.  The  rms  differences  for  stations  1 and  2 (which 
should  be  most  affected  by  waves  in  the  ocean  part  of  the  model)  were 
clearly  reduced  when  waves  were  added  to  the  model. 


119 


LEGEND 


Figure  30.  Prejetty  (1964)  mean  tide  condition  at  gage  1. 
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Figure  31.  Prejetty  (1964)  mean  tide  condition  at  gage  2. 
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Figure  32.  Pre jetty  (1964)  mean  tide  condition  at  gage  3. 
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Figure  34.  Prejetty  (1964)  mean  tide  condition  at  gage  5. 
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Figure  35.  Prejetty  (1964)  mean  tide  condition  at  gage  6. 
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Figure  36.  Prejetty  (1964)  mean  tide  condition  of  tide  in  embayment  area. 
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Figure  42.  Tidal  current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  2S. 
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Figure  44.  Tidal  current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  3C. 
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Figure  45.  Tidal  current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  3S. 


136 


current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  4U. 
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current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  4C. 


Figure  48.  Tidal  current  velocities  for  prejetty  (1964)  and  mean  tide  conditions  at  station  4W 


Tfioco  — o — wro^m 

i i i i i 


poou  qq3 

(s/H)  Aipo|3A 


139 


Figure  49.  Tidal  current  velocities  for  prejetty  (1S64)  and  mean  tide  conditions  at  station  5E. 
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5.  The  Cost  Factor. 


The  solution  of  many  engineering  problems  is  based  not  only  on  the 
reliability  of  the  estimate  but  also  the  cost  and  efforts  involved.  The 
costs  incurred  in  conducting  each  of  the  model  studies  for  Masonboro 
Inlet  are  as  follows: 


The  lumped  parameter  mathematical  model $ 5,000 

(Huval  and  Wintergerst,  1977) 

Mathematical  model 25,400 

(Chen  and  Hembree,  1977) 

Mathematical  model 46,0  35 

(Masch,  Brandes,  and  Reagan,  1977). 

Physical  model $225,000 

(Sager  and  Seabergh,  1977) 


The  cost  for  the  physical  model  represents  an  estimate,  since  the  physi- 
cal model  was  used  to  conduct  other  studies  which  were  not  a part  of  the 
comparison  of  numerical  and  physical  models.  This  shows  that  there  may 
be  a substantial  difference  in  costs  of  carrying  out  tidal  studies  with 
the  various  modeling  techniques.  However,  this  discussion  does  not  cover 
all  aspects  of  the  costs.  Physical  models  are  more  suited  for  some  prob- 
lems than  the  presently  available  numerical  models.  Conversely,  numeri- 
cal models  are  superior  to  physical  models  for  some  purposes.  Numerical 
models  are  newer  than  the  physical  models,  and  the  numerical  technology 
is  advancing  faster.  A judicious  combination  of  physical  and  numerical 
models  will  likely  soon  become  the  most  economical  approach. 

VIII.  SUMMARY,  CONCLUSION,  AND  RECOMMENDATIONS 

1.  Summary. 

To  evaluate  the  state-of-the-art  modeling  techniques  for  tidal  inlets, 
a physical  model,  a lumped  parameter  model,  and  two  two-dimensional  numer- 
ical models  of  Masonboro  Inlet  were  commissioned.  All  four  models  were 
adjusted  to  the  hydrography  of  September  1969,  and  later  used  to  predict 
the  flow  with  prejetty  (1964)  bathymetry.  The  lumped  parameter  model 
could  predict  only  the  average  flow  in  the  inlet  and  the  average  tide 
for  the  basin  at  the  landward  end  of  the  inlet.  The  other  three  models 
yielded  estimates  of  both  tide  and  current  at  several  locations  in  the 
experimental  basin.  All  models  yielded  tide  hydrographs  of  substantially 
the  same  quality.  The  physical  model  and  the  two-dimensional  numerical 
model  developed  by  Masch,  Brandes,  and  Reagan  (1977)  yielded  estimates 
of  currents  of  substantially  equal  quality.  The  lumped  parameter  model 
agreed  reasonably  well  with  the  average  measured  currents  and  inland 
tides  of  1969,  and  with  the  predictions  obtained  with  the  physical  and 
numerical  models  by  Masch,  Brandes,  and  Reagan  for  the  1964  bathymetry. 

The  second  numerical  model  by  Chen  and  Hembree  (1977)  gave  poor  repro- 
duction of  the  measured  currents  of  1969,  and  did  not  agree  well  with 
the  predictions  of  the  other  three  models  for  the  1964  bathymetry. 
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To  extend  the  comparison  of  physical  and  numerical  models,  the  two- 
dimensional  shallow-water  hydrodynamic  equations  were  derived  from  the 
Navier-Stokes  equations  in  Appendix  B.  Terms  describing  the  interaction 
of  horizontal  flows  with  flows  that  are  not  necessarily  horizontal  were 
retained.  The  physical  interpretation  of  each  term  and  the  conditions 
under  which  each  term  can  make  a significant  contribution  to  the  mean 
horizontal  flow  were  discussed  in  Section  II. 

The  design  of  physical  models,  interpreted  as  analog  computers,  and 
the  extent  to  which  the  physical  models  can  reproduce  or  simulate  the 
processes  described  by  the  governing  equations  were  discussed  in  Section 
III.  This  led  to  the  surprising  conclusion  that  the  effects  of  radiation 
stresses,  due  to  wind-generated  waves,  on  the  tidal  flow  can  be  simulated 
with  useful  accuracy  in  a distorted  scale  physical  model.  Sager  and 
Seabergh  (in  preparation,  1977)  independently  show  that  the  agreement 
between  the  Masonboro  Inlet  model  and  the  prototype  could  be  improved 
by  simulating  the  effect  of  waves  on  the  tidal  flow.  This  study  was  not 
a part  of  the  initial  plan  of  investigation  and  the  contribution  of  waves 
was  not  included  in  the  comparison  results  by  any  of  the  models.  It  has 
been  shown  in  this  study  (Sec.  Ill)  that  the  physical  model  could  not 
simulate  the  effects  of  the  earth's  rotation,  the  atmospheric  pressure 
gradient,  or  the  wind  stress;  however,  none  of  these  effects  were  of 
great  importance  in  Masonboro  Inlet.  The  friction  coefficient  in  the 
distorted  scale  Masonboro  Inlet  model  must  be  five  times  as  great  as  the 
bottom  roughness  in  the  prototype  to  correctly  simulate  the  bottom-stress 
term. 

A brief  review  of  the  literature  on  the  numerical  solution  of  the 
two-dimensional  hydrodynamic  equations  indicates  that  the  phase  speed  of 
waves  computed  by  most  numerical  models,  including  the  two  models  used 
for  Masonboro  Inlet,  is  less  than  the  true  speed.  The  ratio  of  the  com- 
puted speed  to  the  true  speed  increases  with  the  ratio  of  wavelength  to 
the  space  increment  (L/AS) . Small-scale  phenomena  with  characteristic 
lengths  less  than  four-space  increments  (4AS)  are  poorly  reproduced  in 
numerical  models.  The  effects  of  the  earth's  rotation,  atmospheric  pres- 
sure gradients,  and  surface  wind  stresses  can  be  adequately  simulated  in 
numerical  models  with  little  difficulty. 

The  two-dimensional  equations  of  motion  were  integrated  in  a cross- 
channel direction  to  form  one-dimensional  channel  equations.  These  in 
turn  were  integrated  along  the  channel  axis  to  form  a zero  dimension  or 
a lumped  parameter  model  of  an  inlet.  This  mode  of  derivation,  where 
the  simplest  equation  is  derived  by  systematically  simplifying  the  most 
complete  equations,  revealed  several  discrepancies  with  the  derivation 
of  the  lumped  parameter  model  discussed  in  this  report. 

In  some  practical  calculations,  these  mathematical  errors  do  not 
lead  to  any  significant  error  in  calculated  tidal  heights  and  currents; 
in  others,  the  errors  may  invalidate  the  technique.  Procedures  for  de- 
termining the  situation  when  the  original  technique  is  valid,  and  for 
extending  the  technique  when  needed,  are  given. 


144 


Both  physical  and  numerical  models  need  to  be  calibrated  for  the 
phenomena  of  interest  because  it  is  impossible  (in  any  one  model)  to 
correctly  model  or  simulate  all  hydrodynamic  phenomena  of  interest  in 
tidal  waterways.  Minor  adjustments  to  even  a well-designed  model  are 
generally  made  during  model  calibration  to  improve  the  agreement  between 
model  and  prototype  for  known  events.  Although  a few  guidelines  have 
been  established  for  the  adjustment  procedure,  adjustment  is  usually  a 
trial-and-error  process. 

2.  Conclusions. 


Both  physical  and  numerical  models  of  inlets,  when  well  designed  and 
properly  used,  can  provide  quantitative  predictions  of  useful  accuracy 
for  hydrodynamic  phenomena.  No  model  can  provide  complete  accuracy  for 
all  phenomena  of  interest.  In  general,  the  physical  model  can  provide 
more  reliable  predictions  than  the  numerical  model  for  the  effect  of 
small-scale  phenomena;  e.g.,  turbulence  and  waves  on  the  phenomena  of 
tidal  scale.  However,  the  numerical  model  can  provide  more  reliable 
predictions  than  physical  models  for  the  effects  of  the  earth's  rotation, 
wind  stress,  and  atmospheric  pressure  gradients. 

The  lumped  parameter  model  developed  by  Keulegan  (1967)  and  many 
extensions  of  this  model,  including  the  model  by  Huval  and  Wintergerst 
(1977),  contain  discrepancies  with  the  rigorous  development  presented  in 
this  report.  This  model  yields  respectable  predictions  for  Masonboro 
Inlet  in  spite  of  the  discrepancies  with  the  rigorous  development  pre- 
sented here,  because  the  processes  incorrectly  described  in  the  model 
are  of  only  minor  importance  in  the  Masonboro  Inlet. 

The  proper  design  or  modification  of  models  and  the  interpretation 
of  results  from  both  physical  and  numerical  models  require  a great  deal 
of  specialized  information  and  training  not  usually  included  in  the  train- 
ing or  experience  of  engineers  or  other  physical  scientists.  This  is 
especially  true  for  numerical  models  which  have  only  recently  become 
common,  and  new  insights  are  gained  almost  every  year.  The  application 
of  models  to  engineering  studies  by  engineers  without  recent  training  or 
experience  in  this  field,  is  unlikely  to  produce  optimal  results  and  the 
results  can  be  misleading. 

When  a numerical  model  capable  of  providing  the  needed  information 
with  only  minor  modification  is  available,  the  model  can  generally  be 
used  more  economically  than  a physical  model  and  can  be  more  easily 
stored  for  future  use. 

3.  Recommendations . 

a.  A systematic  program  should  be  established  for  exploring  .he 
artifacts  of  numerical  modeling  of  tidal  hydraulics  problems,  developing 
solutions  for  such  difficulties  as  the  reduced  speeds  of  short  waves  and 
evaluating  the  relative  merits  of  alternative  procedures.  The  experiences 
of  meteorologists  can  serve  as  guides  for  this  activity,  and  many  of  the 
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results  obtained  can  be  used  with  little  change.  The  boundary  conditions 
needed  for  coastal  models  are  different  from  those  used  in  meteorology; 
therefore,  improved  techniques  are  needed  most  in  the  field  of  boundary 
conditions  and  movable  boundaries. 

b.  Reports  of  modeling  experiments  should  include  candid  discussions 
of  the  specific  efforts  to  improve  model -prototype  calibration  so  that 
data  are  developed  which  will  permit  a systematic  optimization  of  model 
performance. 

c.  Each  modeling  effort  involving  either  physical  or  numerical  models 
should  include  at  least  one  active  team  member  with  previous  training  and 
experience  in  using  the  same  type  of  model. 

d.  Observations  of  wave  direction,  height,  period  or  spectra,  and 
local  winds  should  be  included  with  hydrographic  surveys  in  planning 
data  collection  for  either  physical  or  numerical  models.  The  data 
collection  program  should  include  observations  of  tidal  heights  or  tidal 
currents  at  all  tidal  boundaries  including  any  tidewater  channels  which 
leave  the  model  in  an  inland  direction.  Tidal  observations  for  model 
calibration  and  verification  should  be  continued  for  at  least  1 lunar 
month . 

e.  A lumped  parameter  model  based  on  the  corrected  equations  should 
be  evaluated  for  possible  field  use. 
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APPENDIX  A 


Request  for  Proposals 


TIDAL  CALCULATIONS  FOR  MASONBORO  INLET,  N.  C. 
Introduction 


Masonboro  Inlet,  North  Carolina,  has  been  selected  for  extensive 
field,  physical  model,  and  mathematical  model  investigations  as  a 
"typical"  inlet  problem  area.  As  a basic  part  of  the  investigation,  the 
effectiveness  of  both  mathematic  and  physical  models  in  predicting  tidal 
heights  and  current  velocities  for  selected  prototype  conditions  is  to 
be  determined.  A physical  fixed-bed  model  of  Masonboro  Inlet  has  been 
constructed  and  verified  using  the  data  of  Inclosure  1.  Further  tests 
will  be  conducted  to  arrive  at  the  tidal  heights  and  current  velocities 
in  November  1964  and  June  1967  using  the  inlet  hydrography  and  jetty 
placement  of  those  dates.  Concurrently,  it  is  desired  to  obtain  predic- 
tions using  existing  mathematical  model  procedures  for  the  same  hydro- 
graphic  conditions  to  be  tested  on  the  physical  model.  The  work  of  Earl 
I.  Brown  (1)  and  G.  H.  Keulegan  (2)  is  to  be  considered  as  well  as  other 
appropriate  approaches. 


Purpose 

The  purpose  of  this  study  is  to  evaluate  the  degree  to  which  math- 
ematical models  can  be  used  to  predict  quantitatively  the  hydrodynamics 
(other  than  sediment  movement)  of  flow  through  tidal  inlets.  This  will 
be  done  by  computing  flows  through  a specific  inlet,  Masonboro  Inlet, 
for  which  actual  flow  conditions  are  known. 

The  Prototype 

Masonboro  Inlet  is  located  along  the  North  Carolina  coast  north  of 
the  mouth  of  the  Cape  Fear  River  (see  USC$GS  Chart  No.  1235).  Inclosure 
2 is  the  portion  of  this  chart  showing  the  inlet  area.  The  inlet-bay  is 
an  intricate  complex  of  channels,  overflow  regions,  and  marshland. 

Work  Statement 


The  study  is  to  consist  of  the  calibration  and  application  of  a pre- 
viously selected  lumped  or  discretized  (one- dimensional  or  two- dimens iona 
mathematical  model  to  predict  the  water  surface  time-history  and  current 
velocities  of  Masonboro  Inlet  for  two  specific  hydrographic  conditions. 
Basically,  the  study  will  include:  (a)  adaptation  of  the  previously 
selected  (at  the  proposal  stage)  model  to  Masonboro  Inlet;  (b)  adjustment 
of  the  model  to  allow  reproduction  of  the  prototype  tides  and  currents 
of  12  September  1969  as  shown  by  Inclosure  1;  (c)  prediction  of  the 
tides  and  currents  for  the  additional  hydrographic  conditions  of  the 
inlet  as  shown  by  Inclosures  4 and  5 for  an  idealized  mean  and  spring 


tide  in  the  ocean  (see  Inclosure  6) ; and  (d)  preparation  of  a final 
report  describing  in  detail  all  significant  phases  of  the  study. 


Methodology 


The  proposals  to  be  submitted  in  response  to  the  request  are  to  con- 
sider each  of  the  following  possible  approaches  sepnrately: 

1.  Lumped  approach  (Keulegan) 

2.  Discretized  approach  (one-dimensional) 

3.  Discretized  approach  (two-dimensional) 

Proposals  may  be  submitted  on  any  one,  any  two,  or  all  of  these  approaches 
as  long  as  the  portion  on  each  method  is  kept  completely  separate  so  that 

the  proposal  for  each  separate  method  may  be  accepted  or  rejected  on  its 

own  merits.  Details  of  the  first  method  are  given  in  references  (1)  and 
(2) . The  second  method  may  be  any  existing  model  such  as  that  developed 
by  Reid  and  Bodine  (3)  or  Leendertse  (4).  The  lumped  model  will  be  used 
to  determine  tidal  response  of  the  inner  connecting  channels  and  veloc- 
ities within  the  inlet  due  to  a given  ocean  tide.  The  one-dimensional 
model  will  be  used  to  predict  water  surface  time-histories,  the  tidal 
prism,  and  the  maximum  mean  (average)  velocity  in  the  inlet  channel  and 
interior  channels.  The  two-dimensional  model  will  be  used  to  predict 
variation  of  water  surface  elevations  and  current  velocities  with  time 

in  the  ocean  approach,  the  inlet,  and  the  interior  channels. 


Mathematical  Model  Confirmation 


The  appropriate  topographic  and  field  data  are  attached  as  Inclosure 
1.  Inclosure  3 shows  the  location  of  the  seven  ocean  tide  and  five  cur- 
rent velocity  stations.  The  mathematical  model  will  be  applied  to 
Masonboro  Inlet  using  given  ocean  tide  and  steady  state  or  other  initial 
conditions.  The  model  will  be  adjusted  and  confirmed  to  the  1969  bottom 
survey  and  to  given  1969  velocity  and  tide  data  (Inclosure  1). 

Mathematical  Model  Application 

The  developed  model  will  be  applied  to  Masonboro  Inlet  under  the 
following  conditions: 

1.  Preproject  undeveloped  inlet  conditions — predict  tides,  tidal 
prism,  and  currents  for  November  1964  survey  data  (Inclosure  4)  at 
locations  shown  in  Inclosure  3. 

2.  Modified  inlet  and  north  jetty — predict  tides,  tidal  prism,  and 
currents  for  June  1967  survey  data  (Inclosure  5)  at  locations  shown  in 
Inclosure  3. 
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Report  Preparation 

A final  report  will  be  written  on  the  investigation,  including 
results,  suggestions  for  mathematical  model  usage  appropriate  to  other 
inlets,  and  a description  of  model  development.  It  will  include  a copy 
of  the  program  deck  with  a listing  and  adequate  documentation. 
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Inclosures 

1.  1969  Tidal  Heights,  Current  Velocities,  and  Hydrographic  Survey 

2.  Chart  of  Masonboro  Inlet  Area 

3.  Masonboro  Inlet  Range  and  Gage  Locations 

4.  Hydrographic  Survey  for  November  1964 

5.  Hydrographic  Survey  for  June  1967 

6.  Ocean  Tidal  Elevations 
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MASONBORO  INLET 

PROTOTYPE  VELOCITY  12  SEPT  1969 

RANGE  4 - STATION  W (WEST) 
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Incl  1 to  App.  A,  p.  8 of  15  pp. 
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PROTOTYPE  VELOCITY  12  SEPT  1969 

RANGE  I - STATION  N (NORTH) 
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MEAN  AND  SPRING  TIDES 


APPENDIX  B 


A DERIVATION  OF  THE  LONG  WAVE  EQUATIONS 
1.  General  Considerations. 

Tidal  flows  are  quasi-horizontal  and  nearly  periodic.  Periods  of 
greatest  importance  are  approximately  12  and  24  hours,  but  periods  as 
short  as  3 hours  and  longer  than  24  hours  can  be  identified  in  tidal 
theory  and  in  the  analysis  of  the  tidal  records  from  many  locations. 
Turbulent  eddies,  and  perhaps  wind- generated  waves  with  periods  measured 
in  seconds,  play  a secondary  but  nevertheless  measurable  role  in  modify- 
ing the  flows  which  result  from  astronomical  forces.  Wind- generated 
currents,  though  not  nearly  so  periodic,  may  have  amplitudes  and  dura- 
tions similar  to  the  tides.  Thus,  consideration  of  turbulence,  wind- 
generated waves  and  currents  is  necessary  to  fully  understand  observed 
tidal  flows. 

All  flows  in  the  ocean  or  the  atmosphere  near  sea  level  conform  to 
the  Navier-Stokes  equations  which  express  the  laws  for  the  conservation 
of  momentum  in  a viscous  fluid,  and  to  an  equation  for  the  conservation 
of  mass.  However,  the  Navier-Stokes  equations  are  inconvenient  for  a 
solution  of  the  tidal  hydraulics  problems  because  of  generality.  A con- 
venient set  of  equations  could  be  derived  by  assuming  that  only  the 
motions  of  interest  occur  and  discarding  all  terms  in  the  complete  equa- 
tions that  are  not  obviously  important.  This  procedure  provides  the 
first  approximations  to  a complete  solution,  and  is  adequate  for  many 
applications  but  provides  no  insight  into  the  importance  of  the  neglected 
phenomena.  An  indirect,  but  ultimately  a better  procedure  is  to  filter 
the  undesired  motions  from  the  fundamental  equations  and  retain  the 
terms  which  describe  the  interaction  between  the  motions  of  interest  and 
other  motions  which  may  accompany  them.  The  resulting  equations  may  be 
expressed  in  the  form  of  the  sum  of  the  principal  terms  (which  could 
have  been  derived  directly) , and  groups  of  secondary  terms  which  de- 
scribe the  secondary  phenomena  neglected  in  the  direct  derivation. 

Some  of  the  secondary  phenomena  are  often  highly  correlated  with 
the  primary  flows,  and  effects  on  the  primary  flows  can  often  be  esti- 
mated with  useful  accuracy  from  the  primary  flows.  Hie  nature  of  this 
correlation  can  change  with  changes  in  the  character  of  the  secondary 
flow,  and  may  vanish  altogether.  Identification  of  secondary  phenomena, 
which  may  not  be  modeled  properly,  is  useful  in  comparing  models  for  the 
solution  of  specific  problems  and  for  estimating  the  confidence  which 
should  be  placed  in  a specific  model  for  the  solution  of  a specific 
problem. 

The  equations  of  motion  can  be  simplified  for  astronomical  tides 
and  other  quasi-horizontal  flows  of  similar  scale  by  neglecting  the 
compressibility  of  seawater.  The  incompressibility  approximation  is 
especially  useful  when  numerical  models  are  considered  (see  Sec.  IV). 

The  compressibility  of  seawater  does  not  have  a significant  effect  on 


176 


water  velocities  that  are  less  than  the  speed  of  sound,  or  on  water 
densities  unless  large  vertical  displacements  occur.  Since  tidal  veloc- 
ities are  low  and  the  flow  is  quasi-horizontal,  the  seawater  is  assumed 
incompressible  in  this  derivation. 

Separation  of  tidal  flows  from  small-scale  flows,  such  as  turbulence 
and  wind-generated  waves,  is  more  difficult.  This  separation  is  ac- 
complished, at  least  approximately,  by  defining  a "large-scale  flow" 
which  is  the  average  of  the  total  flow  over  some  element  of  space  and 
time,  and  a "perturbation  flow"  which  is  the  difference  between  the 
total  flow  and  the  large-scale  flow.  Numerical  models  can  only  deal 
effectively  with  the  average  flow  where  large  averaging  volumes  and 
times  are  used  (see  Sec.  IV).  Recognizing  the  difference  in  the 
effective  averaging  implied  by  different  models  is  essential  to  the 
proper  interpretation  of  the  results. 

The  equations  may  be  simplified  by  suppressing  details  of  the  ver- 
tical structure  of  the  flow  by  integration  of  each  term  in  each  equa- 
tion from  the  bottom  to  the  top  of  an  identifiable  fluid  layer. 

The  filtering  techniques  described  here  are  widely  used  in  the  field 
of  hydrodynamics.  The  development  given  below  is  an  outgrowth  of  the 
derivations  by  Fortak  (1962)  and  Harris  (1967). 

2.  The  Primitive  Equations. 

The  equations  for  the  conservation  of  mass  and  momentum  in  an  incom- 
pressible fluid  on  a rotating  earth  may  be  stated  as: 
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where 

v,  w = components  of  velocity  parallel  to  the  x-,  u-,  and  z-axes 
z ■ positive  upward 
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f = Coriolis  parameter  = 2 fi  sin  <f> 

f2  = 2 fl  cos  cj> , a second  Coriolis  parameter 


fi  = angular  velocity  of  the  earth  = 7.292  x 10  5 radians  per  second 
<j>  = latitude 

V = angle  between  the  x-axis  and  the  east  direction 
g = acceleration  of  gravity 
p = fluid  density 

The  components  of  the  stress  tensor  are  given  by: 
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where  a is  the  coefficient  of  molecular  viscosity,  and  p is  the 
pressure.  It  is  necessary  to  consider  the  earth's  rotation  in  these 
equations  because  the  period  of  rotation  is  of  the  same  order  of  magni- 
tude as  the  period  of  the  astronomical  tides.  The  Coriolis  terms  involv- 
ing w are  generally  omitted  because  of  the  quasi-horizontal  character 
of  large-scale  geophysical  flows.  These  terms  are  always  small  relative 
to  the  acceleration  of  gravity,  g,  and  are  omitted  in  this  appendix. 


Treatment  of  the  nonlinear  terms  in  the  momentum  equations  is  simpli- 
fied by  adding  the  product  of  u,  v,  and  w,  and  p times  equation 
(Bl)  to  equation  (B2)  to  obtain: 


u y a o r 

+ U + r— UV  + -r-VM  - fv 

3x  3y  3z 


+ 

r 

,[k  + L. 

l_3t  3x 


3 3 p 

— uv  + — — v^ 
3x  3y 


+ t — vw  + fu 

3 Z 


3 3 2 

uw  + — ^vw  + — -wz  + g 
3y  3z  B 


3 

3x 

raxc  + 

3 

j^y*  + 

3 

3zT2x 

(B5) 

3 

3 

3 

(B6) 

FT1 

lxy  + 

Tpvy  + 

3 

3 

3 

+ 

xz 

3y Xyz 

Sz^sa  • 

(B7) 

3.  Separation  of  Large-  and  Small-Scale  Flows. 


The  unnecessary  details  of  the  small-scale  flow  may  be  eliminated 
from  the  equations  by  partitioning  the  total  flow  into  a mean  value, 
indicated  by  an  overbar  (-) , and  perturbations  around  this  mean  value. 
When  the  upper  boundary  of  the  flow  being  investigated  is  exposed  to  the 
wind,  two  types  of  perturbations  must  be  recognized:  (a)  a generally 
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organized  flow  due  to  the  wind-generated  waves  and  indicated  by  a tilde 
; and  (b)  a generally  random  flow  developed  in  response  to  velocity 
gradients  produced  by  either  the  mean  flow  or  the  waves.  This  random 
perturbation  is  called  turbulence  and  is  indicated  by  a prime  (i).  Thus, 
each  of  the  primary  variables  is  partitioned  into  three  components: 

u = u + u + u'  , v = v + v + v'  , 

w = w + w + w'  , p = p + p + p'  . (B8) 

The  mean  values  are  defined  for  an  elemental  area  or  volume  and  time 
interval  large  enough  to  eliminate  the  variability  due  directly  to  the 
unwanted  small-scale  features  caused  by  waves  and  turbulence,  but  not 
large  enough  to  interfere  with  the  description  of  the  flow  of  primary 
interest. 


Substitution  from  equation  (B8)  into  equations  (Bl)  and  (B5)  to  (B7) 
and  averaging  over  the  same  elemental  space  and  time  interval  will  elim- 
inate linear  terms  in  the  perturbation  quantities,  but  may  leave  products 
of  perturbation  variables.  These  products  of  perturbation  quantities 
reveal  interaction  between  large-  and  small-scale  flows.  The  visible 
wave  on  the  water  surface  is  a result  of  the  velocity  perturbations  due 
to  waves;  thus,  it  is  necessary  to  carry  out  the  vertical  integration 
of  the  equations  before  carrying  out  the  second  averaging  operation  to 
describe  the  complete  interaction  between  the  waves  and  the  mean  flow. 

The  perturbations  resulting  from  turbulence  do  not  have  such  an  identi- 
fiable effect  on  the  elevation  of  the  water  surface  and  the  physics  of 
the  process  is  clearly  revealed  by  completing  the  averaging  process 
for  the  turbulence  terms  before  carrying  out  the  vertical  integration. 
This  procedure  requires  a certain  lack  of  rigor  because  the  character- 
istic periods  of  waves  and  turbulence  overlap,  and  no  completely  rigor- 
ous means  of  distinguishing  between  the  periods  have  been  established. 
This  lack  of  mathematical  rigor  is  believed  to  permit  a more  reliable 
description  of  nature  than  could  be  obtained  by  neglecting  the  most 
prominent  features  of  either  waves  or  turbulence  to  obtain  mathematical 
rigor. 

When  surface  waves  are  omitted,  the  simple  form  of  equations  (Bl) 
and  (B5)  to  (B7)  may  be  recovered  by  combining  the  products  of  turbu- 
lence terms  with  the  definitions  of  the  stresses  given  by  equation  (B4). 
When  treated  in  this  manner,  the  products  of  the  turbulence  terms  are 
generally  called  "Reynolds  stresses",  after  the  British  hydrodynamicist 
who  first  suggested  the  partitioning  process.  The  following  equations 
combine  the  Reynolds  stresses  with  the  molecular  stresses,  and  the  pres- 
sure term  is  expressed  on  the  left  side  of  the  equations,  as  for  inviscid 
fluids.  The  primitive  equations  (Bl)  and  (B6)  to  (B8) , when  modified  by 
the  separation  of  scales  of  motion  and  prepared  for  vertical  integration, 
give: 

|^(u  ♦ u)  + f^(v  ♦ v)  + f^w  + w)  = 0 (B9) 
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where  v = p/p  and  is  assumed  to  be  constant. 
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Each  term  in  equations  (B9)  to  (B12)  must  be  integrated  from  the 
bottom  of  the  fluid  layer,  z = Z(x,y,t),  to  the  top  of  the  layer, 
z = h(x,y,t).  In  a well-mixed  estuary,  the  bottom  of  the  fluid  layer 
is  the  seabed  and  the  top  is  the  free  surface.  In  a stratified  estuary, 
at  least  two  fluid  layers  must  be  considered,  and  the  bottom  of  one  layer 
is  the  top  of  the  next.  Although  this  analysis  is  completed  only  for  a 
single  fluid  layer,  the  general  form  is  used  here  because  it  provides 
symmetry  to  the  development  without  adding  any  complications. 

If  no  fluid  passes  through  the  lower  and  upper  boundaries  of  the 
fluid,  the  boundary  conditions  for  the  vertical  integration  are  given  by: 


3Z  3Z  3Z 
WZ  = 3t  + UZ37  + VZ?y  ’ 


(B14) 
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and 


(BIS) 


where  the  subscript  Z or  h indicates  that  the  variable  is  evaluated 
at  z = Z or  h.  An  overbar  is  used  when  needed  to  show  that  an  average 
value  of  h is  implied.  Z is  assumed  constant  with  respect  to  the 
averaging  process  and  no  fluid  passes  through  the  upper  or  lower 
boundaries. 


When  any  variable,  \p,  is  expressed  in  the  form  ip  = ip  + <p,  the 
average  of  the  vertical  integral  may  be  expressed  as : 

h h h h 

j ip  dz  = J ijT  dz  j (jTdz+  J \p  dz  , (B16) 

Z1  Z1  S Z1 

where  < h.  The  averaging  operator  may  be  taken  inside  the  first  inte- 
gral on  the  right  in  equation  (B16)  because  the  limit  of  integration  is 

constant  with  respect  to  the  averaging  operation.  This  is  not  correct 
for  the  second  and  third  integrals  where  the  upper  limit  has  not  been 
submitted  to  the  averaging  operation.  The  second  overbar  in  the  first 
integral  may  be  omitted  because  ip  is  a constant  with  respect  to  the 
averaging  operation.  Generally,  the  second  overbar  is  not  used. 

The  Leibnitz  rule  for  differentiation  of  an  integral  is  essential  to 
the  formal  integration  of  equations  (B9)  to  (B12) , and  may  be  stated  as: 


ip  dz  - Kh)  H + pm  || . 


(B17) 


Application  of  equations  (B14)  to  (B17)  to  the  continuity  equation 
(B9)  yields:  
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U = j u dz,  V = f v dz  . 

z z 
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U and  V are  vertically  averaged  fluid  transports  parallel  to  the  x- 
and  y-axes. 


The  first  integral  in  equation  (B18)  may  be  expanded  to  form: 


(B20) 


The  average  value  of  the  first  integral  in  equation  (B20)  vanishes 
because  the  upper  limit  is  a constant  with  respect  to  the  averaging  pro- 
cess and  each  term  in  the  integral  is  linear  in  a perturbation  quantity. 
The  second  integral  in  equation  (B20)  may  be  evaluated  by  the  Leibnitz 
rule,  equation  (B17),  to  obtain: 
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(B21) 


The  first  term  in  equation  (B21)  vanishes  because  it  is  the  mean 
value  of  a linear  perturbational  quantity.  The  second  and  third  integrals 
may  be  evaluated  by  expansion  as  a Taylor  series  about  z = h.  Thus,  these 
integrals  take  the  form:  
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where  the  subscript  h indicates  that  the  variables  are  evaluated  at 
z = h.  Note  that  IT  = h - h,  if^j,  v^,  and  their  derivatives  are  also  per- 
turbation quantities,  and  the  mean  value  of  terms  of  odd  order  tends  to 
vanish.  The  first  terms  are  generally  much  larger  than  the  third.  Thus, 
with  good  approximation, equation  (B21)  is  given  by: 
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The  second  integral  in  equation  (B20)  has  a nonzero  value,  even 
though  the  integrand  has  the  form  of  a continuity  equation  which  is  ex- 
pected to  vanish  everywhere  within  the  fluid,  because  the  limit  of  inte- 
gration in  variable  and  is  highly  correlated  with  the  integrand. 

The  third  integral  in  equation  (B18)  may  be  integrated  by  expanding 
each  term  as  a Taylor  series  about  z=  h,  where  h - H = h.  The  average 
value  of  the  resulting  integral  is: 
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All  terms  in  the  above  expression  are  negligibly  small  under  most  con- 
ditions and  are  not  considered  further. 


Equation  (B18)  has  been  derived  on  the  assumption  that  no  fluid 
enters  or  leaves  the  system  through  horizontal  boundaries.  If  fluid  is 

added  at  the  surface  through  precipitation  or  condensation  or  removed 
through  evaporation,  or  if  there  is  seepage  through  the  bottom,  the  net 
effects  of  these  changes  must  be  accounted  for  through  a source-sink 
term.  This  is  indicated  by  R,  on  the  assumption  that  the  maximum  values 
of  the  source-sink  term  will  result  from  rainfall.  Thus,  the  continuity 
equation  may  be  expressed  as: 

!t(S  - z)  4 ih"  * k1  ’ R * E • (B23) 

where 

E = + • (B24) 

E can  be  rewritten  as: 

CB2S) 

where  s is  distance  along  a wave  ray  in  the  direction  of  wave  travel 
and  q^  is  the  velocity  in  the  direction  of  wave  travel  at  z = h,  due 
to  the  wave.  It  results  from  the  correlation  between  the  displacement  of 
free  surface  and  the  horizontal  fluid  velocity  associated  with  the  wave. 
This  transport  contributes  significantly  to  the  horizontal  divergence 
only  in  regions  where  the  wave  amplitude  changes  rapidly  in  the  direction 
of  wave  propagation,  such  as  the  shoaling  and  surf  zones,  or  in  regions, 
such  as  inlets,  where  the  waves  are  modulated  by  the  convergence  or 
divergence  of  the  mean  flow. 


Formal  integration  of  the  equations  for  the  conservation  of  horizon- 
tal momentum  from  the  bottom  to  the  top  of  a fluid  layer  may  be  accom- 
plished through  application  of  equation  (B 14)  to  (B17)  to  obtain: 
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D = h - Z = total  depth. 
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Evaluation  of  equations  (B26)  and  (B27)  requires  the  vertical  inte- 
gral of  the  products  ( 1/p ) 3p/3x  and  (1/p)  3p/3y,  where  p may  be  a 
function  of  z and,  as  yet,  p is  undefined.  An  expression  of  p(z) 
may  be  obtained  by  integrating  the  equation  for  the  conservation  of 
vertical  momentum  from  an  arbitrary  value  of  z to  the  top  of  the  layer. 
The  possible  variation  of  density  with  z may  be  accounted  for  approxi- 
mately by  means  of  the  identity,  ^ 


where 


and 


(B29) 


dz  , 


p"  = p - p . 


Since  a well-mixed  fluid  layer  is  assumed,  it  is  permissible  to  neglect 
(p")2/£3,  and  p"/£>2  may  be  neglected  in  a first  approximation.  Since 
the  error  term  is  retained,  a later  correction  can  possibly  be  estimated 
if  needed. 


In  developing  an  expression  for  the  mean  pressure  through  integration 
of  the  equation  for  vertical  momentum,  it  is  permissible  to  neglect  all 
expressions  involving  the  mean  vertical  velocity  because  of  the  quasi- 
horizontal nature  of  the  mean  flow.  By  combining  equation  (B29)  with 
equations  (B14)  to  (B17)  and  neglecting  the  mean  vertical  velocity, 
equation  (B12)  may  be  integrated  to  obtain: 
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By  rearranging  terms  in  equation  (B30)  and  applying  the  averaging 
operator,  the  desired  expression  for  the  mean  pressure  may  be  obtained 
in  the  form: 


(|)p = + + ez  > 
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(B32) 


The  first  term  on  the  right  in  equation  (B31)  is  the  familiar  hydro- 
static equation  due  to  the  weight  of  the  water  above  z;  the  second 
term  is  the  atmospheric  pressure.  The  third  term  is  a negative  dynamic 
term  due  to  the  vertical  velocity  at  height,  z,  resulting  from  wave 
motion,  and  the  fourth  term,  E2,  represents  closure  terms  which  can  be 
neglected  almost  everywhere. 

It  is  sufficient  to  consider  the  first  integral  in  equation  (B32) 
only  for  the  interval,  h to  h,  because  it  is  a linear  perturbation 
term  and  the  average  value  within  fixed  limits  of  integration  vanishes. 
The  development  of  equation  (B25)  shows  that  the  integral  can  be  obtained 
from  the  first  terms  in  a Taylor  series  expansion  in  the  form  1w.  This 
term  vanishes  almost  everywhere  because  of  the  usual  90°  phase  shift 
between  the  vertical  displacement  of  the  free  surface  and  the  vertical 
velocity  due  to  gravity  waves.  The  mean  values  of  the  second  and  third 
integrals  vanish  almost  everwhere  because  of  the  usual  90°  phase  shift 
between  the  horizontal  and  vertical  velocity  components  due  to  gravity 
waves.  The  fourth  integral  is  always  small  because  of  the  assumption  of 
a well-mixed  fluid  layer. 

In  considering  the  final  integral,  a substitution  for  the  stress 
terms  from  equation  (B13)  will  be  necessary.  When  the  mean  component 
of  the  vertical  velocity  is  neglected, the  required  expression  has  the 
form: 
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z 

The  integrand  in  equation  (B33)  represents  the  gradient  of  the  tur- 
bulent Reynolds  stress  terms  involving  vertical  velocity  and  the 
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horizontal  velocity  in  the  direction  of  flow.  These  Reynolds  stresses 
are  generated  by  the  vertical  shear  in  the  horizontal  mean  flow.  Signif- 
icant values  of  this  gradient  can  occur  only  with  large  values  of  the 
curvature  in  the  mean  flow,  such  as  at  sharp  bends  in  a channel.  The 
term,  (w')^,  represents  the  effect  of  that  part  of  the  roughness  of 
the  free  surface  that  is  generated  by  turbulence  in  the  water.  Signi- 
ficant values  rarely  occur.  The  final  term  is  the  vertical  component 
of  the  turbulent  velocity  at  the  seabed,  and  can  also  be  neglected 
almost  everywhere.  Since  the  pressure  enters  equations  (B26)  and  (B27) 
only  in  the  form  of  horizontal  gradients,  and  all  terms  in  equation 
(B32)  have  been  shown  to  be  small,  equation  (B31)  is  taken  as  an  adequate 
expression  for  p in  the  further  consideration  of  equations  (B26)  and 
(B27).  Substitution  of  the  vertically  averaged  density  p for  p in 
equations  (B26)  and  (B27)  permits  the  application  of  equation  (B16)  to 
be  vertical  integration  of  the  pressure  gradient  terms.  Thus, 
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Z z 
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The  vertical  integral  of  the  pressure  may  be  computed  from  equation  (B31) 
as : 


(?)  I p dz  = (i)  gF  + *2] + (?)  pfcD 

z 

h 

= (j)  gD2  + (j)  W + (t)  V - / w2  dz  , (B36) 

'p'  Z 

where  terms  which  vanish  in  the  mean  have  been  eliminated.  Terms  in  the 
perturbation  of  the  atmospheric  pressure,  p' , due  to  the  weight  of  the 
air  column  between  the  wave  trough  and  the  wave  crest,  are  not  given, 

because  the  ratio  of  air  and  water  densities  is  so  low  as  to  make  the 
perturbation  of  the  atmospheric  pressure  at  the  water  surface  trivial. 


h 

- J w2  dz 

Z 


187 


By  combining  equation  (B36)  with  equations  (B34)  and  (B35) , the  mean 
vertically  integrated  pressure  gradient  terms  are  obtained  in  the  form: 
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(B38) 


In  equations  (B26)  and  (B27) , the  linear  perturbation  terms  contri- 
bute to  the  vertical  integral  only  from  the  region  h to  h,  and  the 
quadratic  terms  contribute  only  from  the  region  Z to  h as  in  the 
derivation  of  equation  (B25) . All  integrals  from  h to  h may  be 
carried  out  through  Taylor  series  expansions  as  with  equation  (B25) . 

By  considering  these  principles  and  using  equations  (B37)  and  (B38) , the 
horizontal  transport  equations  may  be  expressed  as: 


3U 


UV 


— _ fy 

3t  3x  D 3y  P 


+ 


XL J 


XO 


syw  syo 


Sxv  • ^ 

S yv  , CB40) 


where  ( x3X)fj  and  (jzv)h  are  the  mean  horizontal  stresses  parallel  to  the 

° o "o 

x-  and  y-axes  in  the  upper  boundary  of  the  fluid,  (.xzx)h  an^  ^rzy^h  are 

normally  taken  as  the  wind  stress  on  the  water  and  are  evaluated  as 

functions  of  the  wind  velocity,  and  (,tzx)z  and  (jZy) g are  the  mean 

horizontal  stresses  parallel  to  the  x-  and  y-axes  at  the  bottom  of  the 
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fluid  layer.  These  stresses  are  normally  taken  as  the  friction  between 
the  water  and  the  seabed,  and  are  usually  evaluated  as  functions  of  the 
mean  flow  velocity. 

Sxt  end  are  the  molecular  and  turbulent  stresses  experienced  by 
the  interior  of  the  fluid  layer,  where 
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The  molecular  stresses  may  dominate  in  equations  (B41)  and  (B42)  within 
a few  millimeters  of  a solid  surface.  The  turbulent  stresses  are  domi- 
nant in  the  interior  of  the  fluid,  and  act  to  smooth  the  flow  in  a hori- 
zontal plane  in  the  above  equations.  Velocity  shears  due  either  to  the 
mean  flow  or  waves  contribute  to  the  development  of  turbulence. 

sxu  and  syw>  8iven  by: 

h h 

Saw  = 1 7 / (if2  - W2  dz  ♦ ({)  g|^(h) 2 + |^  / uv  dz  , (B43) 

z z 

h _ h_ 
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have  been  called  radiation  stresses  by  Longuet-Higgins  and  Stewart  (1962, 
1963,  1964) , and  are  due  entirely  to  the  Reynolds  stresses  generated  by 
gravity  waves. 

The  terms. 
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are  additional  stresses  resulting  from  waves  that  have  not  been  studied 
in  detail.  Each  of  these  stress  terms  is  proportional  to  the  wave  ampli- 
tude squared.  In  general,  a decrease  in  the  wave  amplitude  with  time  or 
in  the  direction  of  wave  propagation  causes  an  increase  in  the  momentum 
in  the  direction  of  wave  propagation.  These  groups  of  perturbation  terms 
are  generally  smaller  than  the  radiation  stresses,  and  can  generally  be 
neglected  relative  to  the  uncertainties  in  the  radiation  and  turbulence 
stresses.  However,  unrecognized  existence  of  these  terms  could  lead  to 
poor  reproduction  of  radiation  stresses  in  experimental  studies.  The 
final  term  in  each  equation  can  be  neglected  because  of  the  small  pressure 
difference  in  the  atmosphere  between  wave  crest  and  trough,  and  because 
of  the  90°  phase  shift  between  the  displacement  of  the  free  surface  and 
the  horizontal  gradient  of  the  wave-induced  pressure. 

$xv  and  S yv  are  given  by: 
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These  terms  arise  because  the  advective  terms  given  by  equations  (B39) 
and  (B40)  represent  only  the  contribution  from  the  vertically  averaged 
velocity.  In  reality,  the  horizontal  velocity  is  a function  of  z. 

The  current  speed  generally  increases  with  distance  from  the  bottom. 

The  integrals  of  (u")2  and  (v")2  act  as  stresses,  similar  to  the  wind, 
turbulence,  and  radiation  stresses  previously  discussed.  Since  only  the 
gradients  of  these  stresses  affect  the  momentum  balance,  these  terms  can 
play  a significant  role,  only  when  the  gradients  are  large.  In  general, 
these  terms  will  become  large  only  when  the  flow  is  constrained;  e.g., 
when  the  current  passes  through  a channel  connecting  two  bodies  of  water. 

A tidal  inlet  is  such  a channel,  and  these  terms  make  a significant  con- 
tribution to  the  momentum  balance  at  the  entrance  or  exit  of  a tidal 
inlet,  especially  when  the  depth  of  the  inlet  is  less  than  that  of  the 
ocean  or  the  bay  near  the  inlet.  No  measured  data  which  could  be  used  to 
determine  the  importance  of  this  term  are  known. 
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APPENDIX  C 


A DERIVATION  OF  THE  ONE- DIMENSIONAL  OR  CHANNEL  EQUATIONS 


Equations  governing  the  flow  of  quasi-steady  currents  or  long  waves 
in  a channel  may  be  formed  by  choosing  the  x-axis  to  parallel  as  closely 
as  possible  to  the  axis  of  the  channel,  and  to  integrate  each  term  in  the 
long  wave  equations  (1)  and  (3),  with  respect  to  y to  obtain  equations 
governing  the  flow  in  a channel  of  variable  cross  section.  Each  term  in 
the  integrated  equation  must  be  a function  of  x and  t only.  There- 
fore, each  term  is  integrated  or  averaged  over  the  cross  section  for  spe- 
cific values  of  x and  t.  Closure  terms,  generally  neglected  when 
fundamental  equations  are  built  up  from  the  minimum  number  of  terms,  are 
defined.  The  closure  terms  are  negligible  in  many  practical  problems. 
Equations  (1)  and  (B39)  may  be  arranged  and  expressed  as  an  integral  in 
the  form: 
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The  Leibnitz  rule  for  the  differentation  of  an  integral  (eq.  B17) 
can  be  used  to  integrate  the  first  integral  in  equation  (Cl)  in  the 
form: 
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where 


yi 

Ae(x,t)  = J D dy 
yi 

is  the  cross-sectional  area  of  the  inlet  at  position  x along  the  channel 
at  time, t. 


u(x,t) 


y2 

= A31  j U dy  = A31 
yl 


u dz  dy  , 


u"(x,y,z,t)  = u(x,y,z,t)  - u(x,t). 


A double  overbar  (=)  is  used  here  to  indicate  a cross-sectional 
average.  Double  primes  (")  are  used  to  indicate  deviations  from  the 
cross-sectional  average  of  functions  which  vary  with  depth,  or  deviation 
from  the  cross-channel  mean  of  functions  independent  of  z.  Closure 
terms  are  retained  in  compact  form  as,  the  goal  of  this  development  is  to 
show  a generalized,  but  relatively  simple  extension  of  the  lumped  param- 
eter inlet  model  of  Keulegan  (1967). 


The  second  integral  in  equation  (Cl)  may  be  expressed  as: 
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h = b 


1 j h dy,  D = r~~ 


h"  = h - h , D"  = D - D , 


b = b(x)  = y2(x)  - yx(x)  , 


and  yl  (x)  and  y2(x)  are  the  sides  of  the  channel. 
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The  third  integral  in  equation  (Cl)  may  be  expressed  as: 
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where 


(t«c)z  = b_1  J ( ^ zx) Z dy 

yi 

is  the  average  value  of  the  bottom  and  internal  stress  for  an  assigned 
value  of  x and  t. 

The  first  integral  on  the  right  in  equation  (Cl)  may  be  expressed  as: 


2 2 

f j V dy  = fAcv,  v = A'1  J V dy  , 


(C5) 


1 'l 

where  v is  the  average  cross-sectional  value  of  the  flow  normal  to  the 
channel  axis.  This  will  vanish  if  the  cross-channel  flow  is  symmetric 
with  respect  to  the  centerline  of  the  channel. 

The  second  integral  on  the  right  expresses  the  contribution  of  the 
atmosphere  to  the  momentum  of  the  water.  The  variation  of  these  forces 
across  a typical  inlet  for  the  time  periods  of  importance  in  tidal  prob- 
lems is  so  small  that  the  point  values  may  be  treated  as  means,  and  the 
integral  may  be  expressed  as: 


3p 
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The  first  term  on  the  right  in  equation  (C6)  expresses  the  wind 
stress.  Although  wind  stress  cannot  be  predicted  long  in  advance  and  is 
usually  neglected,  it  can  be  significant.  The  second  term,  the  atmos- 
pheric pressure  gradient,  is  nearly  always  negligibly  small. 

The  third  term  on  the  right  in  equation  (Cl)  is  the  cross-channel 
integral  of  the  closure  terms  of  the  two-dimensional  equations  (see 
Sec.  II  and  App.  B) . 

When  v and  the  perturbation  terms  are  neglected,  and  the  relation 
between  bottom  stress  and  velocity  is  considered,  the  only  remaining 
terms  of  potential  significance  in  equation  (2)  and  (B40)  are 


.3h 


= IS)  (Tz^  * * + V - - S; 


y° 
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Division  by  gD  yields: 


3h 
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O/SHVjj-fu  v-s 


ux  ~ ~ 


(C8) 


gD 


gD 


Equation  (C8)  expresses  the  cross-channel  slope  as  a function  01  the 
cross-channel  wind  stress,  the  rotation  of  the  earth,  the  water  depth, 
and  the  closure  terms.  The  equation  is  rarely  used  in  investigations  of 
channel  flow  where  cross-sectional  averages  are  considered.  However, 
it  may  be  useful  in  estimating  the  probable  error  in  using  point  measure- 
ments of  water  surface  elevations  as  estimates  of  the  cross-channel  mean 
surface  elevation. 

The  continuity  equation  for  channel  flow,  obtained  by  cross-channel 
integration  of  equations  (3)  and  (B23) , takes  the  form: 


Ir  ‘ I - / <R  * « ■ 


(C9) 


The  one-dimensional  transport  equation  obtained  by  substitution  from 
equations  (C2)  to  (C6)  into  equation  (Cl)  is: 


ft^3  + jfc"2  + gA4£ + / (T^)z  dy  = Rc 


(CIO) 


where 
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Re  = he  f 6(u"?  + 2uD"U"  + D" (u") 2 dy 
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The  remainder  term,  Rc,  contains  all  terms  which  appear  when  the 
one-dimensional  equations  are  derived,  without  approximation,  from 
the  three-dimensional  equations,  but  which  are  generally  neglected  in 
the  one-dimensional  momentum  equation. 

At  times  it  is  more  convenient  to  work  with  the  equations  for  the 
mean  flow  rather  than  the  equations  for  total  transport.  The  differential 
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equation  for  the  mean  flow  in  a channel  may  be  obtained  from  equation 
(CIO)  by  differentiation  of  the  product  terms  to  obtain: 


. 3u 
^3t  + 


. = 3 3 


. 3h  f 

gA^  + J 


(1sx)2  dy  ♦ 


JMC  3 =1 

UL~  + ^UJ  = 


(C12) 


The  bracketed  term  in  the  above  equation  is  identical  with  the  left 
side  of  equation  (C9).  Thus,  this  term  may  be  eliminated  by  subtracting 
u times  equation  (C9)  from  equation  (C12) . Further  division  by  A yields 
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(R  + E)  dy 
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The  continuity  equation  in  the  proper  form  for  computing  the  mean 
flow  may  be  obtained  by  expressing  KQ  as  Ka  = Db  in  -equation  (C9) , 
differentiating  product  terms  and  dividing  each  term  by  b,  the  width 
of  the  channel,  to  obtain: 


3D  + 3Du 
3t  9x 


-1 
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In  most  applications,  3Z/3t  = 0,  and  the  terms  on  the  right  are 
neglected.  If  also  h <<  Z,  the  equation  may  be  linearized  to  obtain  the 
more  familiar  form: 


3h  3Du  _ 

3t  3x  = 0 ' 


(CIS) 


where  D in  the  case  is  the  difference  between  Z and  the  mean  water 
level . 
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In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 
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Vicksburg,  Miss.  : U.  S.  Waterways  Experiment  Station  ; Spring- 
field,  Va.  : available  from  National  Technical  Information 
Service,  1977. 
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General  investigation  of  tidal  inlets;  a program  of  research 
conducted  jointly  by  U.  S.  Army  Coastal  Engineering  Research 
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